library(knitr)
options(max.print="75")
opts_chunk$set(echo = TRUE, cache = TRUE, warning = FALSE, prompt=FALSE, tidy = TRUE, comment = NA, message = FALSE)
opts_knit$set(width=75)

# Loading library here
library(tidyverse)
#library(ggpubr)
library(MASS)
library(segmented)
library(agricolae)

# set path
source("~/Dropbox/42-JCHWANG-RUTGERS/Projects-Rutgers/Src/utils.R")

Initial setup

Setup color index

# Exposure level
Urban_color <- c("#66c2a5", "#fc8d62", "#8da0cb")
names(Urban_color) <- c("Low", "Medium", "High")
print(Urban_color)
      Low    Medium      High 
"#66c2a5" "#fc8d62" "#8da0cb" 
# plot(1:3, 1:3, col = Urban_color, pch = 16, cex = 4) Ethnicity
Ethnicity_color <- gg_color_hue(3)
names(Ethnicity_color) <- c("SANEMA", "YEKWANA", "Visitors")
print(Ethnicity_color)
   SANEMA   YEKWANA  Visitors 
"#F8766D" "#00BA38" "#619CFF" 
Sanemas <- c("Chajuranha", "Mosenahanha", "Kuyuwininha", "Shianana-Jiyakwanha", "Washudihanha", 
    "Sudukuma")
Yekwana <- c("Kanarakuni", "Fiyakwanha")

Import the data

# Alpha diversities
alpha = readRDS("../data/alpha_imported_all.RDS")
# Mapping This mapping files include all villagers and baseline of visitors of
# all body sites
mt_ms = readRDS("../data/Mapping_MS_Urban.Rds")
Mapping_work <- mt_ms %>% mutate(Urban = ifelse(Ethnicity == "SANEMA", "Low", ifelse(Village == 
    "Fiyakwanha", "Low", ifelse(SampleGroup == "Visitors", "High", "Medium"))) %>% 
    factor(levels = c("Low", "Medium", "High")), Age_num = as.numeric(Age), Age_num = ifelse(is.na(Age_num), 
    as.numeric(gsub(pattern = "_months", replacement = "", Age, ignore.case = T))/12, 
    Age_num), Age_grp1 = ifelse(Age_num >= 18, "Adults", "Children"))
Mapping_work$Age_num[Mapping_work$Age == "4_Days"] <- 0
Mapping_work$Age_num[Mapping_work$Age == "1_day"] <- 0

Set up the Analysis

Work <- merge(alpha, Mapping_work, by = 1) %>% mutate(Age_grp2 = case_when(Age_num <= 
    3 ~ "Age_0-3", Age_num <= 8 ~ "Age_3-8", Age_num < 18 ~ "Age_8-18", T ~ "Adults"))
# Work$Body_Site %>% factor %>% levels
Indices <- colnames(alpha[2:5])
Yrs <- c("2015", "2016")
BodySites <- c("Feces", "Mouth", "Nose", "Right_Arm", "Right_Hand")

Impacts on exposure level

dat_all = list()

Feces

Initial test (no need to run)

We first check the faith_pd which described the overall complexity of the microbial community.

BB = "Feces"
AA = "faith_pd"
# Year 2015
YY = "2015"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers")

# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender + 
    Urban:Gender, data = Work_dat)

# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)

# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num + 
    Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)
Start:  AIC=606.04
get(AA) ~ 1

          Df Sum of Sq    RSS    AIC
+ Urban    1   1050.37 6177.0 583.23
+ Age_num  1    418.01 6809.4 598.63
<none>                 7227.4 606.04
+ Gender   1     88.54 7138.8 606.09

Step:  AIC=583.23
get(AA) ~ Urban

          Df Sum of Sq    RSS    AIC
+ Age_num  1    236.50 5940.5 579.06
+ Gender   1    135.24 6041.8 581.73
<none>                 6177.0 583.23

Step:  AIC=579.06
get(AA) ~ Urban + Age_num

                Df Sum of Sq    RSS    AIC
+ Age_num:Urban  1    177.10 5763.4 576.28
+ Gender         1    173.42 5767.1 576.38
<none>                       5940.5 579.06

Step:  AIC=576.28
get(AA) ~ Urban + Age_num + Urban:Age_num

         Df Sum of Sq    RSS    AIC
+ Gender  1    133.17 5630.2 574.58
<none>                5763.4 576.28

Step:  AIC=574.58
get(AA) ~ Urban + Age_num + Gender + Urban:Age_num

                 Df Sum of Sq    RSS    AIC
<none>                        5630.2 574.58
+ Urban:Gender    1    66.683 5563.5 574.70
+ Age_num:Gender  1    27.428 5602.8 575.81
anova(regforward.out)
Analysis of Variance Table

Response: get(AA)
               Df Sum Sq Mean Sq F value    Pr(>F)    
Urban           1 1050.4 1050.37 28.5436 3.263e-07 ***
Age_num         1  236.5  236.50  6.4268   0.01225 *  
Gender          1  173.4  173.42  4.7126   0.03148 *  
Urban:Age_num   1  136.9  136.86  3.7190   0.05565 .  
Residuals     153 5630.2   36.80                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The best model from 2015 includes Urban, Age, Urban:Age interaction and Gender (in decreasing importance).

# Year 2016
YY = "2016"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers")

# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender + 
    Urban:Gender, data = Work_dat)

# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)

# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num + 
    Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)
Start:  AIC=655.24
get(AA) ~ 1

          Df Sum of Sq    RSS    AIC
+ Urban    1    703.20 5886.2 636.70
+ Age_num  1    219.36 6370.0 651.07
<none>                 6589.4 655.24
+ Gender   1     10.27 6579.1 656.95

Step:  AIC=636.7
get(AA) ~ Urban

          Df Sum of Sq    RSS    AIC
+ Age_num  1   136.141 5750.1 634.44
<none>                 5886.2 636.70
+ Gender   1    27.951 5858.3 637.83

Step:  AIC=634.44
get(AA) ~ Urban + Age_num

                Df Sum of Sq    RSS    AIC
+ Age_num:Urban  1   170.585 5579.5 630.96
<none>                       5750.1 634.44
+ Gender         1    37.584 5712.5 635.24

Step:  AIC=630.96
get(AA) ~ Urban + Age_num + Urban:Age_num

         Df Sum of Sq    RSS    AIC
<none>                5579.5 630.96
+ Gender  1    20.061 5559.4 632.30
anova(regforward.out)
Analysis of Variance Table

Response: get(AA)
               Df Sum Sq Mean Sq F value    Pr(>F)    
Urban           1  703.2  703.20 22.4340 4.428e-06 ***
Age_num         1  136.1  136.14  4.3433   0.03858 *  
Urban:Age_num   1  170.6  170.59  5.4421   0.02077 *  
Residuals     178 5579.5   31.35                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The best model from 2016 includes Urban, Age, Urban:Age interaction (in decreasing importance), but Gender was not there.

So we will fitting the data with Urban, Age, Urban:Age interaction, and Gender for both years to keep unity.

# Final model
lm_models <- list()
for (YY in Yrs) {
    Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers")
    fit.final <- lm(as.formula(paste0(AA, " ~ Age_num + Urban + Age_num:Urban + Gender")), 
        data = Work_dat)
    print(paste0("--------Below is ANOVA table for year ", YY))
    anova(fit.final) %>% print()
    
    print(paste0("--------Below is coefficient for year ", YY))
    summary(fit.final) %>% print()
    lm_models[[YY]] <- fit.final
}
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: faith_pd
               Df Sum Sq Mean Sq F value    Pr(>F)    
Age_num         1  418.0  418.01 11.3594 0.0009503 ***
Urban           1  868.9  868.86 23.6110 2.893e-06 ***
Gender          1  173.4  173.42  4.7126 0.0314847 *  
Age_num:Urban   1  136.9  136.86  3.7190 0.0556480 .  
Residuals     153 5630.2   36.80                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = as.formula(paste0(AA, " ~ Age_num + Urban + Age_num:Urban + Gender")), 
    data = Work_dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-14.1948  -4.6478   0.4276   4.8227  14.0085 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)         26.79492    1.04178  25.720  < 2e-16 ***
Age_num              0.04259    0.03749   1.136   0.2578    
UrbanMedium         -7.04666    1.47745  -4.769 4.27e-06 ***
GenderM             -1.87206    0.98408  -1.902   0.0590 .  
Age_num:UrbanMedium  0.12776    0.06625   1.928   0.0556 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.066 on 153 degrees of freedom
Multiple R-squared:  0.221, Adjusted R-squared:  0.2006 
F-statistic: 10.85 on 4 and 153 DF,  p-value: 9.041e-08

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: faith_pd
               Df Sum Sq Mean Sq F value    Pr(>F)    
Age_num         1  219.4  219.36  6.9839  0.008962 ** 
Urban           1  620.0  619.99 19.7390 1.561e-05 ***
Gender          1   37.6   37.58  1.1966  0.275488    
Age_num:Urban   1  153.1  153.06  4.8732  0.028563 *  
Residuals     177 5559.4   31.41                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = as.formula(paste0(AA, " ~ Age_num + Urban + Age_num:Urban + Gender")), 
    data = Work_dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-16.7396  -3.4067  -0.0396   3.6630  16.7508 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)         24.03558    0.84472  28.454  < 2e-16 ***
Age_num              0.02299    0.02872   0.800   0.4246    
UrbanMedium         -6.36322    1.38225  -4.604  7.9e-06 ***
GenderM             -0.68035    0.85130  -0.799   0.4253    
Age_num:UrbanMedium  0.13204    0.05982   2.208   0.0286 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.604 on 177 degrees of freedom
Multiple R-squared:  0.1563,    Adjusted R-squared:  0.1372 
F-statistic: 8.198 on 4 and 177 DF,  p-value: 4.349e-06

Next we perform the broken stick regression based on the selected final model

for (YY in Yrs) {
    # YY = '2016'
    fit.seg <- segmented(lm_models[[YY]], seg.Z = ~Age_num, psi = c(3))
    summary(fit.seg)
    plot(fit.seg$fitted.values, fit.seg$residuals)
    plot(fit.seg$model$faith_pd, fit.seg$fitted.values)
    plot(fit.seg$model$Age_num, hatvalues(fit.seg))
    plot(fit.seg$model$Age_num, cooks.distance(fit.seg))
    plot(hatvalues(fit.seg), fit.seg$residuals)
}

It seems from the residual plots that the point with faith_pd smaller than 10 is dragging the curve down and causing non-random residual.

Work %>% filter(Body_Site == BB, faith_pd < 12, SampleGroup == "Villagers")
  Row.names faith_pd  pielou_e  shannon observed_otus BarcodeSequence
1   X15.104 11.52192 0.6547718 3.985899            68    ATTTAGGACGAC
2   X15.144  8.21380 0.5725987 3.248024            51    GGTTCCATTAGG
3   X15.149 10.49707 0.6060992 3.455032            52    CGATAGGCCTTA
    LinkerPrimerSequence Barcode_location Tube_ID Year SampleGroup Body_Site
1 CCGGACTACHVGGGTWTCTAAT          P03_02C  15-104 2015   Villagers     Feces
2 CCGGACTACHVGGGTWTCTAAT          P03_03D  15-144 2015   Villagers     Feces
3 CCGGACTACHVGGGTWTCTAAT          P03_03F  15-149 2015   Villagers     Feces
  Body_Site_Type    Village Ethnicity                Date    House_ID Family_ID
1          Fecal Kanarakuni   YEKWANA 10/1/2015-10/3/2015     2015_13    2015_4
2          Fecal Kanarakuni   YEKWANA 10/1/2015-10/3/2015 2015_38_39c   2015_21
3          Fecal Kanarakuni   YEKWANA 10/1/2015-10/3/2015    2015_6_7   2015_17
  Subject_ID Gender      Age  Urban Age_num Age_grp1 Age_grp2
1    2015_25      M        1 Medium    1.00 Children  Age_0-3
2   2015_108      M       15 Medium   15.00 Children Age_8-18
3   2015_143      M 9_months Medium    0.75 Children  Age_0-3
 [ reached 'max' / getOption("max.print") -- omitted 10 rows ]

It does appears that many baby are outliers, so re-do the analysis without baby.

Modeling without babies

dat_p = list()
dat_p[["2015"]] = list()
dat_p[["2016"]] = list()
faith_pd

Overview

BB = "Feces"
AA = "faith_pd"

Work_dat <- Work %>% filter(Body_Site == BB, SampleGroup == "Villagers", Age_num >= 
    1)
p <- ggplot(Work_dat, aes(x = Age_num, y = get(AA), color = Urban, shape = Gender)) + 
    geom_point(size = 1.5) + scale_shape_manual(values = c(16, 21)) + scale_color_manual(breaks = names(Urban_color), 
    values = Urban_color) + facet_grid(. ~ Year) + theme_jw() + theme(aspect.ratio = 0.8, 
    panel.background = element_rect(color = "black")) + labs(x = "Age (Years)", y = AA, 
    title = BB)
print(p)

Modeling 2015

# Year 2015
YY = "2015"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
    Age_num >= 1)

# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender + 
    Urban:Gender, data = Work_dat)

# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)

# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num + 
    Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)

anova(regforward.out)

Modeling 2016

# Year 2016
YY = "2016"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
    Age_num >= 1)

# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender + 
    Urban:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)

# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num + 
    Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)

anova(regforward.out)

Final linear model

# Final linear model
lm_models <- list()
for (YY in Yrs) {
    Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
        Age_num >= 1)
    fit.final <- lm(paste0(AA, " ~ Age_num + Urban + Gender"), data = Work_dat)
    print(paste0("--------Below is ANOVA table for year ", YY))
    anova(fit.final) %>% print()
    
    print(paste0("--------Below is coefficient for year ", YY))
    summary(fit.final) %>% print()
    lm_models[[YY]] <- fit.final
}
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: faith_pd
           Df Sum Sq Mean Sq F value    Pr(>F)    
Age_num     1  251.8  251.81  7.2084  0.008073 ** 
Urban       1  760.7  760.71 21.7765 6.735e-06 ***
Gender      1  207.1  207.05  5.9272  0.016083 *  
Residuals 150 5239.9   34.93                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Age_num + Urban + Gender"), data = Work_dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-15.2172  -4.9553   0.3526   4.6727  12.2899 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 26.70408    0.96592  27.646  < 2e-16 ***
Age_num      0.06583    0.03033   2.170   0.0316 *  
UrbanMedium -4.71461    0.98345  -4.794  3.9e-06 ***
GenderM     -2.34055    0.96138  -2.435   0.0161 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.91 on 150 degrees of freedom
Multiple R-squared:  0.1888,    Adjusted R-squared:  0.1726 
F-statistic: 11.64 on 3 and 150 DF,  p-value: 6.705e-07

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: faith_pd
           Df Sum Sq Mean Sq F value    Pr(>F)    
Age_num     1   74.2   74.17  2.6078 0.1081605    
Urban       1  436.0  436.01 15.3293 0.0001297 ***
Gender      1   62.2   62.18  2.1861 0.1410845    
Residuals 173 4920.7   28.44                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Age_num + Urban + Gender"), data = Work_dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-14.4679  -3.4498  -0.1077   3.4365  16.8511 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 24.21478    0.78318  30.919  < 2e-16 ***
Age_num      0.03226    0.02421   1.333    0.184    
UrbanMedium -3.49149    0.87185  -4.005 9.21e-05 ***
GenderM     -1.20086    0.81220  -1.479    0.141    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.333 on 173 degrees of freedom
Multiple R-squared:  0.1042,    Adjusted R-squared:  0.08866 
F-statistic: 6.708 on 3 and 173 DF,  p-value: 0.0002618

Next we perform the broken stick regression based on the selected final model

stick_models <- list()
for (YY in Yrs) {
    # YY = '2016'
    Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
        Age_num >= 1)
    fit.seg <- segmented(lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
    summary(fit.seg) %>% print()
    plot(fit.seg$fitted.values, fit.seg$residuals)
    plot(fit.seg$model$faith_pd, fit.seg$fitted.values)
    plot(fit.seg$model$Age_num, hatvalues(fit.seg))
    plot(fit.seg$model$Age_num, cooks.distance(fit.seg))
    plot(hatvalues(fit.seg), fit.seg$residuals)
    anova(fit.seg, lm_models[[YY]]) %>% print()
    davies.test(lm_models[[YY]], seg.Z = ~Age_num, k = 10) %>% print()
    confint.segmented(fit.seg) %>% print()
    qqnorm(fit.seg$residuals)
    slope(fit.seg) %>% print()
    stick_models[[YY]] <- fit.seg
}

    ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)

Estimated Break-Point(s):
               Est. St.Err
psi1.Age_num 7.918   1.69

Meaningful coefficients of the linear terms:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  20.2457     2.0376   9.936  < 2e-16 ***
Age_num       1.1480     0.4147   2.768  0.00636 ** 
UrbanMedium  -4.3660     0.9399  -4.645 7.45e-06 ***
GenderM      -2.2691     0.9107  -2.492  0.01382 *  
U1.Age_num   -1.1658     0.4157  -2.804       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.578 on 148 degrees of freedom
Multiple R-Squared: 0.287,  Adjusted R-squared: 0.2629 

Convergence attained in 2 iter. (rel. change 0)

Analysis of Variance Table

Model 1: faith_pd ~ Age_num + Urban + Gender + U1.Age_num + psi1.Age_num
Model 2: faith_pd ~ Age_num + Urban + Gender
  Res.Df    RSS Df Sum of Sq     F    Pr(>F)    
1    148 4605.7                                 
2    150 5239.9 -2   -634.19 10.19 7.146e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    Davies' test for a change in the slope

data:  formula = faith_pd ~ Age_num + Urban + Gender ,   method = lm 
model = gaussian , link = identity  
segmented variable = Age_num
'best' at = 7.5556, n.points = 9, p-value = 0.0001387
alternative hypothesis: two.sided

                Est. CI(95%).low CI(95%).up
psi1.Age_num 7.91839     4.57902    11.2578

$Age_num
            Est.  St.Err.  t value CI(95%).l CI(95%).u
slope1  1.148000 0.414700  2.76820  0.328470  1.967500
slope2 -0.017832 0.035963 -0.49584 -0.088899  0.053236


    ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)

Estimated Break-Point(s):
              Est. St.Err
psi1.Age_num   28  5.359

Meaningful coefficients of the linear terms:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 22.55375    0.94600  23.841  < 2e-16 ***
Age_num      0.18250    0.06320   2.888  0.00438 ** 
UrbanMedium -3.80783    0.85894  -4.433 1.66e-05 ***
GenderM     -1.09063    0.79021  -1.380  0.16933    
U1.Age_num  -0.31478    0.09054  -3.477       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.184 on 171 degrees of freedom
Multiple R-Squared: 0.1634,  Adjusted R-squared: 0.1389 

Convergence attained in 3 iter. (rel. change 8.2738e-07)

Analysis of Variance Table

Model 1: faith_pd ~ Age_num + Urban + Gender + U1.Age_num + psi1.Age_num
Model 2: faith_pd ~ Age_num + Urban + Gender
  Res.Df    RSS Df Sum of Sq      F   Pr(>F)   
1    171 4595.5                                
2    173 4920.7 -2   -325.12 6.0488 0.002896 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    Davies' test for a change in the slope

data:  formula = faith_pd ~ Age_num + Urban + Gender ,   method = lm 
model = gaussian , link = identity  
segmented variable = Age_num
'best' at = 29.889, n.points = 9, p-value = 0.005568
alternative hypothesis: two.sided

                Est. CI(95%).low CI(95%).up
psi1.Age_num 27.9995     17.4216    38.5775

$Age_num
           Est.  St.Err. t value CI(95%).l  CI(95%).u
slope1  0.18250 0.063195  2.8879   0.05776  0.3072500
slope2 -0.13228 0.064494 -2.0510  -0.25959 -0.0049704

Creating dataset for plots

dat15 <- Work %>% filter(Body_Site == BB, Year == "2015", SampleGroup == "Villagers", 
    Age_num >= 1)
dat16 <- Work %>% filter(Body_Site == BB, Year == "2016", SampleGroup == "Villagers", 
    Age_num >= 1)

dat15 <- cbind(dat15, predict(stick_models[["2015"]], newdata = dat15, interval = "confidence"))
dat16 <- cbind(dat16, predict(stick_models[["2016"]], newdata = dat16, interval = "confidence"))

dat_p[["2015"]][[AA]] = dat15
dat_p[["2016"]][[AA]] = dat16
observed_otu

Overview

BB = "Feces"
AA = "observed_otus"

Work_dat <- Work %>% filter(Body_Site == BB, SampleGroup == "Villagers", Age_num >= 
    1)
p <- ggplot(Work_dat, aes(x = Age_num, y = get(AA), color = Urban, shape = Gender)) + 
    geom_point(size = 1.5) + scale_shape_manual(values = c(16, 21)) + scale_color_manual(breaks = names(Urban_color), 
    values = Urban_color) + facet_grid(. ~ Year) + theme_jw() + theme(aspect.ratio = 0.8, 
    panel.background = element_rect(color = "black")) + labs(x = "Age (Years)", y = AA, 
    title = BB)
print(p)

# Year 2015
YY = "2015"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
    Age_num >= 1)

# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender + 
    Urban:Gender, data = Work_dat)

# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)

# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num + 
    Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)

anova(regforward.out)
# Year 2016
YY = "2016"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
    Age_num >= 1)

# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender + 
    Urban:Gender, data = Work_dat)

# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)

# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num + 
    Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)

anova(regforward.out)
# Final linear model
lm_models <- list()
for (YY in Yrs) {
    Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
        Age_num >= 1)
    fit.final <- lm(paste0(AA, " ~ Age_num + Urban + Gender"), data = Work_dat)
    print(paste0("--------Below is ANOVA table for year ", YY))
    anova(fit.final) %>% print()
    
    print(paste0("--------Below is coefficient for year ", YY))
    summary(fit.final) %>% print()
    lm_models[[YY]] <- fit.final
}
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: observed_otus
           Df Sum Sq Mean Sq F value    Pr(>F)    
Age_num     1  41607   41607  6.7746 0.0101742 *  
Urban       1  92047   92047 14.9875 0.0001611 ***
Gender      1  33787   33787  5.5013 0.0203122 *  
Residuals 150 921238    6142                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Age_num + Urban + Gender"), data = Work_dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-166.815  -60.028   -7.596   56.206  200.962 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 241.7114    12.8076  18.872  < 2e-16 ***
Age_num       0.8844     0.4022   2.199 0.029421 *  
UrbanMedium -52.0959    13.0400  -3.995 0.000101 ***
GenderM     -29.8986    12.7473  -2.345 0.020312 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 78.37 on 150 degrees of freedom
Multiple R-squared:  0.1538,    Adjusted R-squared:  0.1369 
F-statistic: 9.088 on 3 and 150 DF,  p-value: 1.449e-05

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: observed_otus
           Df Sum Sq Mean Sq F value   Pr(>F)   
Age_num     1   3124    3124  0.8962 0.345123   
Urban       1  32964   32964  9.4556 0.002447 **
Gender      1   3173    3173  0.9100 0.341438   
Residuals 173 603117    3486                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Age_num + Urban + Gender"), data = Work_dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-155.109  -46.839   -0.734   43.746  195.489 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 205.3456     8.6706  23.683  < 2e-16 ***
Age_num       0.1914     0.2680   0.714  0.47600    
UrbanMedium -30.2261     9.6523  -3.131  0.00204 ** 
GenderM      -8.5779     8.9919  -0.954  0.34144    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 59.04 on 173 degrees of freedom
Multiple R-squared:  0.06112,   Adjusted R-squared:  0.04484 
F-statistic: 3.754 on 3 and 173 DF,  p-value: 0.01204

Next we perform the broken stick regression based on the selected final model

stick_models <- list()
for (YY in Yrs) {
    # YY = '2015'
    Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
        Age_num >= 1)
    fit.seg <- segmented(lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
    summary(fit.seg) %>% print()
    plot(fit.seg$fitted.values, fit.seg$residuals)
    plot(fit.seg$model$observed_otus, fit.seg$fitted.values)
    anova(fit.seg, lm_models[[YY]]) %>% print()
    davies.test(lm_models[[YY]], seg.Z = ~Age_num, k = 10) %>% print()
    confint.segmented(fit.seg) %>% print()
    qqnorm(fit.seg$residuals)
    slope(fit.seg) %>% print()
    stick_models[[YY]] <- fit.seg
}

    ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)

Estimated Break-Point(s):
               Est. St.Err
psi1.Age_num 8.333  1.747

Meaningful coefficients of the linear terms:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  170.512     24.597   6.932 1.19e-10 ***
Age_num       12.481      4.284   2.913 0.004131 ** 
UrbanMedium  -48.451     12.635  -3.835 0.000186 ***
GenderM      -29.088     12.246  -2.375 0.018818 *  
U1.Age_num   -12.605      4.307  -2.927       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 75.27 on 148 degrees of freedom
Multiple R-Squared: 0.2298,  Adjusted R-squared: 0.2038 

Convergence attained in 2 iter. (rel. change 0)

Analysis of Variance Table

Model 1: observed_otus ~ Age_num + Urban + Gender + U1.Age_num + psi1.Age_num
Model 2: observed_otus ~ Age_num + Urban + Gender
  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1    148 838484                                  
2    150 921238 -2    -82755 7.3035 0.0009442 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    Davies' test for a change in the slope

data:  formula = observed_otus ~ Age_num + Urban + Gender ,   method = lm 
model = gaussian , link = identity  
segmented variable = Age_num
'best' at = 7.5556, n.points = 9, p-value = 0.001613
alternative hypothesis: two.sided

               Est. CI(95%).low CI(95%).up
psi1.Age_num 8.3333     4.88077    11.7858

$Age_num
           Est. St.Err.  t value CI(95%).l CI(95%).u
slope1 12.48100 4.28410  2.91330    4.0150  20.94700
slope2 -0.12465 0.52149 -0.23902   -1.1552   0.90588


    ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)

Estimated Break-Point(s):
                Est. St.Err
psi1.Age_num 28.004  6.004

Meaningful coefficients of the linear terms:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 188.2708    10.4042  18.096  < 2e-16 ***
Age_num       1.7358     0.6676   2.600 0.010131 *  
UrbanMedium -33.4771     9.6220  -3.479 0.000638 ***
GenderM      -7.4450     8.8148  -0.845 0.399509    
U1.Age_num   -3.2363     1.0119  -3.198       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 57.67 on 171 degrees of freedom
Multiple R-Squared: 0.1146,  Adjusted R-squared: 0.08873 

Convergence attained in 5 iter. (rel. change 4.3811e-06)

Analysis of Variance Table

Model 1: observed_otus ~ Age_num + Urban + Gender + U1.Age_num + psi1.Age_num
Model 2: observed_otus ~ Age_num + Urban + Gender
  Res.Df    RSS Df Sum of Sq     F   Pr(>F)   
1    171 568752                               
2    173 603117 -2    -34365 5.166 0.006631 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    Davies' test for a change in the slope

data:  formula = observed_otus ~ Age_num + Urban + Gender ,   method = lm 
model = gaussian , link = identity  
segmented variable = Age_num
'best' at = 29.889, n.points = 9, p-value = 0.01317
alternative hypothesis: two.sided

                Est. CI(95%).low CI(95%).up
psi1.Age_num 28.0042     16.1526    39.8559

$Age_num
          Est. St.Err. t value CI(95%).l CI(95%).u
slope1  1.7358 0.66756  2.6002   0.41809    3.0535
slope2 -1.5005 0.76088 -1.9721  -3.00250    0.0014

Creating data for plots

dat15 <- Work %>% filter(Body_Site == BB, Year == "2015", SampleGroup == "Villagers", 
    Age_num >= 1)
dat16 <- Work %>% filter(Body_Site == BB, Year == "2016", SampleGroup == "Villagers", 
    Age_num >= 1)

dat15 <- cbind(dat15, predict(stick_models[["2015"]], newdata = dat15, interval = "confidence"))
dat16 <- cbind(dat16, predict(stick_models[["2016"]], newdata = dat16, interval = "confidence"))

# lm model children
for (YY in Yrs) {
    Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
        Age_num >= 1, Age_num <= 8)
    fit.child <- lm(paste0(AA, " ~ Age_num + Urban + Gender"), data = Work_dat)
    print(paste0("--------Below is ANOVA table for year ", YY))
    anova(fit.child) %>% print()
    
    print(paste0("--------Below is coefficient for year ", YY))
    summary(fit.child) %>% print()
}
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: observed_otus
          Df Sum Sq Mean Sq F value    Pr(>F)    
Age_num    1  68041   68041 14.3091 0.0003902 ***
Urban      1  10433   10433  2.1940 0.1443576    
Gender     1   9950    9950  2.0924 0.1538140    
Residuals 54 256776    4755                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Age_num + Urban + Gender"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-111.68  -47.16  -15.23   32.11  198.99 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  154.237     25.038   6.160 9.43e-08 ***
Age_num       13.557      3.991   3.396  0.00129 ** 
UrbanMedium  -26.009     18.732  -1.388  0.17069    
GenderM      -26.505     18.324  -1.447  0.15381    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 68.96 on 54 degrees of freedom
Multiple R-squared:  0.2562,    Adjusted R-squared:  0.2148 
F-statistic: 6.198 on 3 and 54 DF,  p-value: 0.00107

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: observed_otus
          Df Sum Sq Mean Sq F value   Pr(>F)   
Age_num    1   3546  3546.3  1.2597 0.266855   
Urban      1  22170 22169.6  7.8753 0.007036 **
Gender     1    489   488.6  0.1736 0.678674   
Residuals 52 146385  2815.1                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Age_num + Urban + Gender"), data = Work_dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-108.472  -37.501   -5.733   28.342  110.281 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  181.711     20.126   9.028 3.14e-12 ***
Age_num        3.761      3.575   1.052  0.29760    
UrbanMedium  -43.276     15.427  -2.805  0.00706 ** 
GenderM       -6.009     14.424  -0.417  0.67867    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 53.06 on 52 degrees of freedom
Multiple R-squared:  0.1518,    Adjusted R-squared:  0.1029 
F-statistic: 3.103 on 3 and 52 DF,  p-value: 0.03442
dat_p[["2015"]][[AA]] = dat15
dat_p[["2016"]][[AA]] = dat16
shannon
BB = "Feces"
AA = "shannon"

Work_dat <- Work %>% filter(Body_Site == BB, SampleGroup == "Villagers", Age_num >= 
    1)
p <- ggplot(Work_dat, aes(x = Age_num, y = get(AA), color = Urban, shape = Gender)) + 
    geom_point(size = 1.5) + scale_shape_manual(values = c(16, 21)) + scale_color_manual(breaks = names(Urban_color), 
    values = Urban_color) + facet_grid(. ~ Year) + theme_jw() + theme(aspect.ratio = 0.8, 
    panel.background = element_rect(color = "black")) + labs(x = "Age (Years)", y = AA, 
    title = BB)
print(p)

# Year 2015
YY = "2015"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
    Age_num >= 1)

# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender + 
    Urban:Gender, data = Work_dat)

# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)

# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num + 
    Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)

anova(regforward.out)
# Year 2016
YY = "2016"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
    Age_num >= 1)

# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender + 
    Urban:Gender, data = Work_dat)

# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)

# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num + 
    Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)

anova(regforward.out)
# Final linear model
lm_models <- list()
for (YY in Yrs) {
    Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
        Age_num >= 1)
    fit.final <- lm(paste0(AA, " ~ Age_num + Urban + Gender"), data = Work_dat)
    print(paste0("--------Below is ANOVA table for year ", YY))
    anova(fit.final) %>% print()
    
    print(paste0("--------Below is coefficient for year ", YY))
    summary(fit.final) %>% print()
    lm_models[[YY]] <- fit.final
}
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: shannon
           Df Sum Sq Mean Sq F value   Pr(>F)   
Age_num     1  0.778  0.7776  1.1775 0.279607   
Urban       1  5.371  5.3707  8.1321 0.004962 **
Gender      1  2.583  2.5833  3.9115 0.049788 * 
Residuals 150 99.064  0.6604                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Age_num + Urban + Gender"), data = Work_dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.18307 -0.55013  0.05821  0.61292  1.88644 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  5.822744   0.132813  43.842  < 2e-16 ***
Age_num      0.003427   0.004171   0.822  0.41256    
UrbanMedium -0.399803   0.135222  -2.957  0.00361 ** 
GenderM     -0.261435   0.132188  -1.978  0.04979 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.8127 on 150 degrees of freedom
Multiple R-squared:  0.081, Adjusted R-squared:  0.06262 
F-statistic: 4.407 on 3 and 150 DF,  p-value: 0.005302

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: shannon
           Df Sum Sq Mean Sq F value    Pr(>F)    
Age_num     1  0.304  0.3042  0.5727 0.4502035    
Urban       1  7.123  7.1228 13.4092 0.0003325 ***
Gender      1  0.185  0.1845  0.3473 0.5563921    
Residuals 173 91.896  0.5312                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Age_num + Urban + Gender"), data = Work_dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.23807 -0.46918 -0.03367  0.54295  1.72129 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  5.552453   0.107028  51.879  < 2e-16 ***
Age_num      0.001486   0.003308   0.449 0.653878    
UrbanMedium -0.439992   0.119146  -3.693 0.000297 ***
GenderM     -0.065415   0.110994  -0.589 0.556392    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7288 on 173 degrees of freedom
Multiple R-squared:  0.07649,   Adjusted R-squared:  0.06048 
F-statistic: 4.776 on 3 and 173 DF,  p-value: 0.003181

Next we perform the broken stick regression based on the selected final model

stick_models <- list()
for (YY in Yrs) {
    # YY = '2015'
    Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
        Age_num >= 1)
    fit.seg <- segmented(lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
    summary(fit.seg) %>% print()
    plot(fit.seg$fitted.values, fit.seg$residuals)
    plot(fit.seg$model$shannon, fit.seg$fitted.values)
    anova(fit.seg, lm_models[[YY]]) %>% print()
    davies.test(lm_models[[YY]], seg.Z = ~Age_num, k = 10) %>% print()
    confint.segmented(fit.seg) %>% print()
    qqnorm(fit.seg$residuals)
    slope(fit.seg) %>% print()
    stick_models[[YY]] <- fit.seg
}

    ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)

Estimated Break-Point(s):
               Est. St.Err
psi1.Age_num 3.571  1.078

Meaningful coefficients of the linear terms:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   4.7605     0.4957   9.605   <2e-16 ***
Age_num       0.3248     0.2050   1.584   0.1152    
UrbanMedium  -0.3437     0.1352  -2.542   0.0120 *  
GenderM      -0.2380     0.1309  -1.818   0.0711 .  
U1.Age_num   -0.3252     0.2050  -1.586       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7994 on 148 degrees of freedom
Multiple R-Squared: 0.1227,  Adjusted R-squared: 0.09308 

Convergence attained in 2 iter. (rel. change 4.5034e-16)

Analysis of Variance Table

Model 1: shannon ~ Age_num + Urban + Gender + U1.Age_num + psi1.Age_num
Model 2: shannon ~ Age_num + Urban + Gender
  Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
1    148 94.567                              
2    150 99.064 -2   -4.4971 3.5191 0.03213 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    Davies' test for a change in the slope

data:  formula = shannon ~ Age_num + Urban + Gender ,   method = lm 
model = gaussian , link = identity  
segmented variable = Age_num
'best' at = 7.5556, n.points = 9, p-value = 0.1001
alternative hypothesis: two.sided

               Est. CI(95%).low CI(95%).up
psi1.Age_num 3.5712     1.44061     5.7018

$Age_num
              Est.   St.Err.   t value CI(95%).l CI(95%).u
slope1  0.32478000 0.2049900  1.584400 -0.080303 0.7298700
slope2 -0.00043922 0.0044672 -0.098321 -0.009267 0.0083885


    ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)

Estimated Break-Point(s):
                Est. St.Err
psi1.Age_num 25.999 11.682

Meaningful coefficients of the linear terms:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  5.441128   0.139278  39.067  < 2e-16 ***
Age_num      0.011852   0.009988   1.187 0.237032    
UrbanMedium -0.461718   0.120135  -3.843 0.000171 ***
GenderM     -0.057888   0.111406  -0.520 0.604003    
U1.Age_num  -0.019809   0.012968  -1.528       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.728 on 171 degrees of freedom
Multiple R-Squared: 0.08932,  Adjusted R-squared: 0.06269 

Convergence *not* attained in 1 iter. (rel. change 9.0543e-07)

Analysis of Variance Table

Model 1: shannon ~ Age_num + Urban + Gender + U1.Age_num + psi1.Age_num
Model 2: shannon ~ Age_num + Urban + Gender
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    171 90.619                           
2    173 91.896 -2   -1.2765 1.2044 0.3024

    Davies' test for a change in the slope

data:  formula = shannon ~ Age_num + Urban + Gender ,   method = lm 
model = gaussian , link = identity  
segmented variable = Age_num
'best' at = 22.667, n.points = 9, p-value = 0.5631
alternative hypothesis: two.sided

                Est. CI(95%).low CI(95%).up
psi1.Age_num 25.9988     2.94015    49.0575

$Age_num
             Est.   St.Err.  t value  CI(95%).l CI(95%).u
slope1  0.0118520 0.0099883  1.18660 -0.0078642 0.0315680
slope2 -0.0079572 0.0081803 -0.97273 -0.0241050 0.0081901

Create data for plots

dat15 <- Work %>% filter(Body_Site == BB, Year == "2015", SampleGroup == "Villagers", 
    Age_num >= 1)
dat16 <- Work %>% filter(Body_Site == BB, Year == "2016", SampleGroup == "Villagers", 
    Age_num >= 1)

dat15 <- cbind(dat15, predict(stick_models[["2015"]], newdata = dat15, interval = "confidence"))
dat16 <- cbind(dat16, predict(stick_models[["2016"]], newdata = dat16, interval = "confidence"))


# lm model children
for (YY in Yrs) {
    Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
        Age_num >= 1, Age_num <= 8)
    fit.child <- lm(paste0(AA, " ~ Age_num + Urban + Gender"), data = Work_dat)
    print(paste0("--------Below is ANOVA table for year ", YY))
    anova(fit.child) %>% print()
    
    print(paste0("--------Below is coefficient for year ", YY))
    summary(fit.child) %>% print()
}
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: shannon
          Df Sum Sq Mean Sq F value   Pr(>F)   
Age_num    1  4.928  4.9281  8.1245 0.006172 **
Urban      1  0.487  0.4867  0.8023 0.374378   
Gender     1  0.094  0.0945  0.1558 0.694633   
Residuals 54 32.755  0.6066                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Age_num + Urban + Gender"), data = Work_dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.79544 -0.55073  0.02717  0.42599  1.58893 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  5.06029    0.28279  17.894   <2e-16 ***
Age_num      0.11652    0.04508   2.585   0.0125 *  
UrbanMedium -0.18395    0.21157  -0.869   0.3884    
GenderM     -0.08168    0.20695  -0.395   0.6946    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7788 on 54 degrees of freedom
Multiple R-squared:  0.144, Adjusted R-squared:  0.09642 
F-statistic: 3.028 on 3 and 54 DF,  p-value: 0.03722

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: shannon
          Df  Sum Sq Mean Sq F value  Pr(>F)  
Age_num    1  0.2542 0.25419  0.5781 0.45050  
Urban      1  3.0389 3.03886  6.9109 0.01124 *
Gender     1  0.0061 0.00609  0.0138 0.90677  
Residuals 52 22.8653 0.43972                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Age_num + Urban + Gender"), data = Work_dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.87772 -0.38368 -0.04977  0.32432  1.43925 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  5.65450    0.25154  22.479   <2e-16 ***
Age_num     -0.03775    0.04468  -0.845   0.4020    
UrbanMedium -0.50693    0.19281  -2.629   0.0112 *  
GenderM      0.02122    0.18028   0.118   0.9068    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6631 on 52 degrees of freedom
Multiple R-squared:  0.1261,    Adjusted R-squared:  0.07567 
F-statistic: 2.501 on 3 and 52 DF,  p-value: 0.06955
dat_p[["2015"]][[AA]] = dat15
dat_p[["2016"]][[AA]] = dat16
Data for feces
dat_all[[BB]] = dat_p

Mouth

dat_p = list()
dat_p[["2015"]] = list()
dat_p[["2016"]] = list()

Modeling without babies

faith_pd

Overview

BB = "Mouth"
AA = "faith_pd"

Work_dat <- Work %>% filter(Body_Site == BB, SampleGroup == "Villagers", Age_num >= 
    1)
p <- ggplot(Work_dat, aes(x = Age_num, y = get(AA), color = Urban, shape = Gender)) + 
    geom_point(size = 1.5) + scale_shape_manual(values = c(16, 21)) + scale_color_manual(breaks = names(Urban_color), 
    values = Urban_color) + facet_grid(. ~ Year) + theme_jw() + theme(aspect.ratio = 0.8, 
    panel.background = element_rect(color = "black")) + labs(x = "Age (Years)", y = AA, 
    title = BB)
print(p)

# Year 2015
YY = "2015"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
    Age_num >= 1)

# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender + 
    Urban:Gender, data = Work_dat)

# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)

# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num + 
    Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)

anova(regforward.out)

For year 2015, only urban and age is significant in determing the faith_pd

# Year 2016
YY = "2016"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
    Age_num >= 1)

# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender + 
    Urban:Gender, data = Work_dat)

# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)

# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num + 
    Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)

anova(regforward.out)

In 2016, In addition to Age_num and Urban, the interaction term between age_num and urban is also significant at 0.05. For consistency, only age_num and urban were included.

# Final linear model
lm_models <- list()
for (YY in Yrs) {
    Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
        Age_num >= 1)
    fit.final <- lm(paste0(AA, " ~ Age_num + Urban"), data = Work_dat)
    print(paste0("--------Below is ANOVA table for year ", YY))
    anova(fit.final) %>% print()
    
    print(paste0("--------Below is coefficient for year ", YY))
    summary(fit.final) %>% print()
    lm_models[[YY]] <- fit.final
}
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: faith_pd
           Df  Sum Sq Mean Sq F value    Pr(>F)    
Age_num     1  326.20  326.20  17.846 3.998e-05 ***
Urban       1  219.03  219.03  11.983 0.0006877 ***
Residuals 161 2942.88   18.28                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Age_num + Urban"), data = Work_dat)

Residuals:
   Min     1Q Median     3Q    Max 
-6.638 -2.297 -0.327  1.292 35.558 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 12.01104    0.63355  18.958  < 2e-16 ***
Age_num      0.07478    0.02165   3.454 0.000707 ***
UrbanMedium -2.38385    0.68866  -3.462 0.000688 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.275 on 161 degrees of freedom
Multiple R-squared:  0.1563,    Adjusted R-squared:  0.1458 
F-statistic: 14.91 on 2 and 161 DF,  p-value: 1.142e-06

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: faith_pd
           Df  Sum Sq Mean Sq F value    Pr(>F)    
Age_num     1  578.15  578.15  35.856 1.616e-08 ***
Urban       1  172.54  172.54  10.700   0.00134 ** 
Residuals 144 2321.87   16.12                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Age_num + Urban"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.7790 -2.8766  0.1378  2.7790 12.5137 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 11.07857    0.59478  18.626  < 2e-16 ***
Age_num      0.11555    0.02022   5.714 6.12e-08 ***
UrbanMedium -2.24572    0.68652  -3.271  0.00134 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.015 on 144 degrees of freedom
Multiple R-squared:  0.2443,    Adjusted R-squared:  0.2338 
F-statistic: 23.28 on 2 and 144 DF,  p-value: 1.739e-09

Next we perform the broken stick regression based on the selected final model

stick_models <- list()
for (YY in Yrs) {
    # YY = '2015'
    Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
        Age_num >= 1)
    fit.seg <- segmented(lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
    summary(fit.seg) %>% print()
    plot(fit.seg$fitted.values, fit.seg$residuals)
    plot(fit.seg$model$faith_pd, fit.seg$fitted.values)
    qqnorm(fit.seg$residuals)
    
    if (!"segmented" %in% class(fit.seg)) {
        print(paste(YY, BB, AA, paste(YY, BB, AA, "same as lm model", sep = "-"), 
            sep = "-"))
    } else {
        anova(fit.seg, lm_models[[YY]]) %>% print()
        davies.test(lm_models[[YY]], seg.Z = ~Age_num, k = 10) %>% print()
        confint.segmented(fit.seg) %>% print()
        stick_models[[YY]] <- fit.seg
        slope(fit.seg) %>% print()
    }
    
    
}

    ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)

Estimated Break-Point(s):
                Est. St.Err
psi1.Age_num 17.874  9.239

Meaningful coefficients of the linear terms:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 12.89173    0.92518  13.934  < 2e-16 ***
Age_num     -0.02790    0.09509  -0.293 0.769591    
UrbanMedium -2.38755    0.68826  -3.469 0.000672 ***
U1.Age_num   0.15067    0.10410   1.447       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.268 on 159 degrees of freedom
Multiple R-Squared: 0.1697,  Adjusted R-squared: 0.1488 

Convergence attained in 5 iter. (rel. change 1.57e-16)

Analysis of Variance Table

Model 1: faith_pd ~ Age_num + Urban + U1.Age_num + psi1.Age_num
Model 2: faith_pd ~ Age_num + Urban
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    159 2896.3                           
2    161 2942.9 -2   -46.535 1.2773 0.2816

    Davies' test for a change in the slope

data:  formula = faith_pd ~ Age_num + Urban ,   method = lm 
model = gaussian , link = identity  
segmented variable = Age_num
'best' at = 20.667, n.points = 9, p-value = 0.3251
alternative hypothesis: two.sided

                Est. CI(95%).low CI(95%).up
psi1.Age_num 17.8742   -0.371911    36.1202
$Age_num
           Est.  St.Err.  t value CI(95%).l CI(95%).u
slope1 -0.02790 0.095089 -0.29341 -0.215700   0.15990
slope2  0.12277 0.043318  2.83430  0.037222   0.20833


    ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)

Estimated Break-Point(s):
                Est. St.Err
psi1.Age_num 64.428  6.366

Meaningful coefficients of the linear terms:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 10.95981    0.60095  18.237  < 2e-16 ***
Age_num      0.12518    0.02154   5.811 3.92e-08 ***
UrbanMedium -2.30255    0.68742  -3.350  0.00104 ** 
U1.Age_num  -0.82889    0.94576  -0.876       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.011 on 142 degrees of freedom
Multiple R-Squared: 0.2563,  Adjusted R-squared: 0.2354 

Convergence *not* attained in 36 iter. (rel. change -0.004225)

Analysis of Variance Table

Model 1: faith_pd ~ Age_num + Urban + U1.Age_num + psi1.Age_num
Model 2: faith_pd ~ Age_num + Urban
  Res.Df    RSS Df Sum of Sq     F Pr(>F)
1    142 2285.1                          
2    144 2321.9 -2   -36.817 1.144 0.3215

    Davies' test for a change in the slope

data:  formula = faith_pd ~ Age_num + Urban ,   method = lm 
model = gaussian , link = identity  
segmented variable = Age_num
'best' at = 66, n.points = 9, p-value = 0.3474
alternative hypothesis: two.sided

                Est. CI(95%).low CI(95%).up
psi1.Age_num 64.4275      51.844     77.011
$Age_num
           Est.  St.Err.  t value CI(95%).l CI(95%).u
slope1  0.12518 0.021542  5.81090  0.082596   0.16777
slope2 -0.70370 0.945510 -0.74426 -2.572800   1.16540

For mouth, broken stick regression did not significantly better than a simple multiple linear regression. So following graph is made using lm.

Data for plots

dat15 <- Work %>% filter(Body_Site == BB, Year == "2015", SampleGroup == "Villagers", 
    Age_num >= 1)
dat16 <- Work %>% filter(Body_Site == BB, Year == "2016", SampleGroup == "Villagers", 
    Age_num >= 1)

dat15 <- cbind(dat15, predict(lm_models[["2015"]], newdata = dat15, interval = "confidence"))
dat16 <- cbind(dat16, predict(lm_models[["2016"]], newdata = dat16, interval = "confidence"))

dat_p[["2015"]][[AA]] = dat15
dat_p[["2016"]][[AA]] = dat16
observed_otu

Overview

BB = "Mouth"
AA = "observed_otus"

Work_dat <- Work %>% filter(Body_Site == BB, SampleGroup == "Villagers", Age_num >= 
    1)
p <- ggplot(Work_dat, aes(x = Age_num, y = get(AA), color = Urban, shape = Gender)) + 
    geom_point(size = 1.5) + scale_shape_manual(values = c(16, 21)) + scale_color_manual(breaks = names(Urban_color), 
    values = Urban_color) + facet_grid(. ~ Year) + theme_jw() + theme(aspect.ratio = 0.8, 
    panel.background = element_rect(color = "black")) + labs(x = "Age (Years)", y = AA, 
    title = BB)
print(p)

# Year 2015
YY = "2015"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
    Age_num >= 1)

# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender + 
    Urban:Gender, data = Work_dat)

# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)

# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num + 
    Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)

anova(regforward.out)
# Year 2016
YY = "2016"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
    Age_num >= 1)

# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender + 
    Urban:Gender, data = Work_dat)

# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)

# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num + 
    Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)

anova(regforward.out)
# Final linear model
lm_models <- list()
for (YY in Yrs) {
    Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
        Age_num >= 1)
    fit.final <- lm(paste0(AA, " ~ Age_num + Urban"), data = Work_dat)
    print(paste0("--------Below is ANOVA table for year ", YY))
    anova(fit.final) %>% print()
    
    print(paste0("--------Below is coefficient for year ", YY))
    summary(fit.final) %>% print()
    lm_models[[YY]] <- fit.final
}
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: observed_otus
           Df Sum Sq Mean Sq F value    Pr(>F)    
Age_num     1   9449  9449.1  7.0084 0.0089190 ** 
Urban       1  20971 20970.6 15.5540 0.0001194 ***
Residuals 161 217068  1348.2                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Age_num + Urban"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-55.840 -24.481  -4.688  14.206 210.453 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  87.1007     5.4412  16.008  < 2e-16 ***
Age_num       0.3370     0.1859   1.812 0.071811 .  
UrbanMedium -23.3257     5.9144  -3.944 0.000119 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 36.72 on 161 degrees of freedom
Multiple R-squared:  0.1229,    Adjusted R-squared:  0.112 
F-statistic: 11.28 on 2 and 161 DF,  p-value: 2.6e-05

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: observed_otus
           Df Sum Sq Mean Sq F value    Pr(>F)    
Age_num     1  28639 28638.6 18.8704 2.625e-05 ***
Urban       1  10545 10545.4  6.9485  0.009309 ** 
Residuals 144 218541  1517.6                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Age_num + Urban"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-77.785 -28.554   0.636  27.055 114.412 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  78.1266     5.7704  13.539  < 2e-16 ***
Age_num       0.8092     0.1962   4.125 6.25e-05 ***
UrbanMedium -17.5569     6.6604  -2.636  0.00931 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 38.96 on 144 degrees of freedom
Multiple R-squared:  0.152, Adjusted R-squared:  0.1403 
F-statistic: 12.91 on 2 and 144 DF,  p-value: 6.968e-06

Next we perform the broken stick regression based on the selected final model

stick_models <- list()
for (YY in Yrs) {
    # YY = '2015'
    Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
        Age_num >= 1)
    fit.seg <- segmented(lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
    summary(fit.seg) %>% print()
    plot(fit.seg$fitted.values, fit.seg$residuals)
    plot(fit.seg$model$observed_otus, fit.seg$fitted.values)
    qqnorm(fit.seg$residuals)
    
    if (!"segmented" %in% class(fit.seg)) {
        print(paste(YY, BB, AA, "same as lm model", sep = "-"))
    } else {
        anova(fit.seg, lm_models[[YY]]) %>% print()
        davies.test(lm_models[[YY]], seg.Z = ~Age_num, k = 10) %>% print()
        confint.segmented(fit.seg) %>% print()
        stick_models[[YY]] <- fit.seg
        slope(fit.seg) %>% print()
    }
}

    ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)

Estimated Break-Point(s):
                Est. St.Err
psi1.Age_num 23.714  8.224

Meaningful coefficients of the linear terms:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  94.6545     6.9218  13.675  < 2e-16 ***
Age_num      -0.4631     0.5465  -0.847 0.398024    
UrbanMedium -23.0323     5.8751  -3.920 0.000131 ***
U1.Age_num    1.5114     0.7270   2.079       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 36.45 on 159 degrees of freedom
Multiple R-Squared: 0.1463,  Adjusted R-squared: 0.1248 

Convergence attained in 1 iter. (rel. change 8.3498e-07)

Analysis of Variance Table

Model 1: observed_otus ~ Age_num + Urban + U1.Age_num + psi1.Age_num
Model 2: observed_otus ~ Age_num + Urban
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    159 211292                           
2    161 217068 -2   -5775.9 2.1732 0.1172

    Davies' test for a change in the slope

data:  formula = observed_otus ~ Age_num + Urban ,   method = lm 
model = gaussian , link = identity  
segmented variable = Age_num
'best' at = 14.111, n.points = 9, p-value = 0.1755
alternative hypothesis: two.sided

                Est. CI(95%).low CI(95%).up
psi1.Age_num 23.7139     7.47148    39.9563
$Age_num
           Est. St.Err.  t value CI(95%).l CI(95%).u
slope1 -0.46312 0.54649 -0.84744 -1.542400    0.6162
slope2  1.04830 0.48321  2.16940  0.093923    2.0026


    ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)

Estimated Break-Point(s):
              Est. St.Err
psi1.Age_num   20 11.171

Meaningful coefficients of the linear terms:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  85.60221    8.99102   9.521   <2e-16 ***
Age_num      -0.02298    0.93448  -0.025   0.9804    
UrbanMedium -15.89449    6.74892  -2.355   0.0199 *  
U1.Age_num    1.26284    1.01494   1.244       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 38.93 on 142 degrees of freedom
Multiple R-Squared: 0.1651,  Adjusted R-squared: 0.1416 

Convergence attained in 1 iter. (rel. change 4.0401e-10)

Analysis of Variance Table

Model 1: observed_otus ~ Age_num + Urban + U1.Age_num + psi1.Age_num
Model 2: observed_otus ~ Age_num + Urban
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    142 215168                           
2    144 218541 -2   -3373.4 1.1132 0.3314

    Davies' test for a change in the slope

data:  formula = observed_otus ~ Age_num + Urban ,   method = lm 
model = gaussian , link = identity  
segmented variable = Age_num
'best' at = 22.667, n.points = 9, p-value = 0.7082
alternative hypothesis: two.sided

             Est. CI(95%).low CI(95%).up
psi1.Age_num   20    -2.08335    42.0833
$Age_num
            Est. St.Err.   t value CI(95%).l CI(95%).u
slope1 -0.022982 0.93448 -0.024593  -1.87030    1.8243
slope2  1.239900 0.38531  3.217800   0.47816    2.0015

Again, linear model is sufficiently good

Data for plots

dat15 <- Work %>% filter(Body_Site == BB, Year == "2015", SampleGroup == "Villagers", 
    Age_num >= 1)
dat16 <- Work %>% filter(Body_Site == BB, Year == "2016", SampleGroup == "Villagers", 
    Age_num >= 1)

dat15 <- cbind(dat15, predict(lm_models[["2015"]], newdata = dat15, interval = "confidence"))
dat16 <- cbind(dat16, predict(lm_models[["2016"]], newdata = dat16, interval = "confidence"))

dat_p[["2015"]][[AA]] = dat15
dat_p[["2016"]][[AA]] = dat16
shannon
BB = "Mouth"
AA = "shannon"

Work_dat <- Work %>% filter(Body_Site == BB, SampleGroup == "Villagers", Age_num >= 
    1)
p <- ggplot(Work_dat, aes(x = Age_num, y = get(AA), color = Urban, shape = Gender)) + 
    geom_point(size = 1.5) + scale_shape_manual(values = c(16, 21)) + scale_color_manual(breaks = names(Urban_color), 
    values = Urban_color) + facet_grid(. ~ Year) + theme_jw() + theme(aspect.ratio = 0.8, 
    panel.background = element_rect(color = "black")) + labs(x = "Age (Years)", y = AA, 
    title = BB)
print(p)

# Year 2015
YY = "2015"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
    Age_num >= 1)

# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender + 
    Urban:Gender, data = Work_dat)

# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)

# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num + 
    Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)

anova(regforward.out)
# Year 2016
YY = "2016"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
    Age_num >= 1)

# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender + 
    Urban:Gender, data = Work_dat)

# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)

# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num + 
    Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)

anova(regforward.out)
# Final linear model
lm_models <- list()
for (YY in Yrs) {
    Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
        Age_num >= 1)
    fit.final <- lm(paste0(AA, " ~ Age_num + Urban"), data = Work_dat)
    print(paste0("--------Below is ANOVA table for year ", YY))
    anova(fit.final) %>% print()
    
    print(paste0("--------Below is coefficient for year ", YY))
    summary(fit.final) %>% print()
    lm_models[[YY]] <- fit.final
}
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: shannon
           Df  Sum Sq Mean Sq F value    Pr(>F)    
Age_num     1   1.100  1.1004  1.5312    0.2177    
Urban       1  29.613 29.6126 41.2071 1.466e-09 ***
Residuals 161 115.699  0.7186                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Age_num + Urban"), data = Work_dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.71659 -0.53916  0.02071  0.46772  3.11725 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.4484667  0.1256208  35.412  < 2e-16 ***
Age_num     -0.0002609  0.0042929  -0.061    0.952    
UrbanMedium -0.8765321  0.1365469  -6.419 1.47e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.8477 on 161 degrees of freedom
Multiple R-squared:  0.2098,    Adjusted R-squared:    0.2 
F-statistic: 21.37 on 2 and 161 DF,  p-value: 5.877e-09

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: shannon
           Df  Sum Sq Mean Sq F value  Pr(>F)  
Age_num     1   3.651  3.6509  4.0000 0.04738 *
Urban       1   0.109  0.1088  0.1192 0.73044  
Residuals 144 131.430  0.9127                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Age_num + Urban"), data = Work_dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.52186 -0.54464  0.02496  0.70104  1.85779 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.782972   0.141509  26.733   <2e-16 ***
Age_num      0.009464   0.004812   1.967   0.0511 .  
UrbanMedium -0.056385   0.163336  -0.345   0.7304    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9554 on 144 degrees of freedom
Multiple R-squared:  0.02781,   Adjusted R-squared:  0.01431 
F-statistic:  2.06 on 2 and 144 DF,  p-value: 0.1312

Next we perform the broken stick regression based on the selected final model

stick_models <- list()
for (YY in Yrs) {
    # YY = '2015'
    Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
        Age_num >= 1)
    fit.seg <- segmented(lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
    summary(fit.seg) %>% print()
    plot(fit.seg$fitted.values, fit.seg$residuals)
    plot(fit.seg$model$shannon, fit.seg$fitted.values)
    qqnorm(fit.seg$residuals)
    
    if (!"segmented" %in% class(fit.seg)) {
        print(paste(YY, BB, AA, "same as lm model", sep = "-"))
    } else {
        anova(fit.seg, lm_models[[YY]]) %>% print()
        davies.test(lm_models[[YY]], seg.Z = ~Age_num, k = 10) %>% print()
        confint.segmented(fit.seg) %>% print()
        stick_models[[YY]] <- fit.seg
        slope(fit.seg) %>% print()
    }
    
    
}

Call:
lm(formula = paste0(AA, " ~ Age_num + Urban"), data = Work_dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.71659 -0.53916  0.02071  0.46772  3.11725 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.4484667  0.1256208  35.412  < 2e-16 ***
Age_num     -0.0002609  0.0042929  -0.061    0.952    
UrbanMedium -0.8765321  0.1365469  -6.419 1.47e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.8477 on 161 degrees of freedom
Multiple R-squared:  0.2098,    Adjusted R-squared:    0.2 
F-statistic: 21.37 on 2 and 161 DF,  p-value: 5.877e-09

[1] "2015-Mouth-shannon-same as lm model"

    ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)

Estimated Break-Point(s):
               Est. St.Err
psi1.Age_num 3.867 11.041

Meaningful coefficients of the linear terms:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.98951    0.91249   4.372 2.36e-05 ***
Age_num     -0.04860    0.37329  -0.130    0.897    
UrbanMedium -0.04946    0.16576  -0.298    0.766    
U1.Age_num   0.05857    0.37336   0.157       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9617 on 142 degrees of freedom
Multiple R-Squared: 0.02855,  Adjusted R-squared: 0.00119 

Convergence attained in 3 iter. (rel. change -2.2215e-07)

Analysis of Variance Table

Model 1: shannon ~ Age_num + Urban + U1.Age_num + psi1.Age_num
Model 2: shannon ~ Age_num + Urban
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    142 131.33                           
2    144 131.43 -2  -0.10071 0.0544  0.947

    Davies' test for a change in the slope

data:  formula = shannon ~ Age_num + Urban ,   method = lm 
model = gaussian , link = identity  
segmented variable = Age_num
'best' at = 22.667, n.points = 9, p-value = 0.7103
alternative hypothesis: two.sided

                Est. CI(95%).low CI(95%).up
psi1.Age_num 3.86673    -17.9592    25.6927
$Age_num
             Est.   St.Err. t value   CI(95%).l CI(95%).u
slope1 -0.0486030 0.3732900 -0.1302 -0.78653000  0.689320
slope2  0.0099667 0.0051258  1.9444 -0.00016605  0.020099

Again broken stick is not sufficiently better than linear

data for plots

dat15 <- Work %>% filter(Body_Site == BB, Year == "2015", SampleGroup == "Villagers", 
    Age_num >= 1)
dat16 <- Work %>% filter(Body_Site == BB, Year == "2016", SampleGroup == "Villagers", 
    Age_num >= 1)

dat15 <- cbind(dat15, predict(lm_models[["2015"]], newdata = dat15, interval = "confidence"))
dat16 <- cbind(dat16, predict(lm_models[["2016"]], newdata = dat16, interval = "confidence"))

dat_p[["2015"]][[AA]] = dat15
dat_p[["2016"]][[AA]] = dat16
Data for mouth
dat_all[[BB]] = dat_p

Nose

dat_p = list()
dat_p[["2015"]] = list()
dat_p[["2016"]] = list()

Modeling without babies

faith_pd
BB = "Nose"
AA = "faith_pd"

Work_dat <- Work %>% filter(Body_Site == BB, SampleGroup == "Villagers", Age_num >= 
    1)
p <- ggplot(Work_dat, aes(x = Age_num, y = get(AA), color = Urban, shape = Gender)) + 
    geom_point(size = 1.5) + scale_shape_manual(values = c(16, 21)) + scale_color_manual(breaks = names(Urban_color), 
    values = Urban_color) + facet_grid(. ~ Year) + theme_jw() + theme(aspect.ratio = 0.8, 
    panel.background = element_rect(color = "black")) + labs(x = "Age (Years)", y = AA, 
    title = BB)
print(p)

# Year 2015
YY = "2015"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
    Age_num >= 1)

# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender + 
    Urban:Gender, data = Work_dat)

# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)

# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num + 
    Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)

anova(regforward.out)

For year 2015, only urban and age is significant in determing the faith_pd

# Year 2016
YY = "2016"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
    Age_num >= 1)

# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender + 
    Urban:Gender, data = Work_dat)

# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)

# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num + 
    Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)

anova(regforward.out)

In 2016, Age_num and Urban as well.

For consistency, only age_num and urban were included.

# Final linear model
lm_models <- list()
for (YY in Yrs) {
    Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
        Age_num >= 1)
    fit.final <- lm(paste0(AA, " ~ Age_num + Urban"), data = Work_dat)
    print(paste0("--------Below is ANOVA table for year ", YY))
    anova(fit.final) %>% print()
    
    print(paste0("--------Below is coefficient for year ", YY))
    summary(fit.final) %>% print()
    lm_models[[YY]] <- fit.final
}
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: faith_pd
           Df  Sum Sq Mean Sq F value   Pr(>F)   
Age_num     1   784.6  784.55  4.5739 0.034276 * 
Urban       1  1613.6 1613.59  9.4072 0.002615 **
Residuals 134 22984.6  171.53                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Age_num + Urban"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-18.705  -9.132  -3.084   7.123  44.949 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 13.64549    2.04108   6.685 5.67e-10 ***
Age_num      0.18238    0.06913   2.638  0.00932 ** 
UrbanMedium  7.18087    2.34124   3.067  0.00262 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 13.1 on 134 degrees of freedom
Multiple R-squared:  0.09448,   Adjusted R-squared:  0.08096 
F-statistic: 6.991 on 2 and 134 DF,  p-value: 0.001295

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: faith_pd
           Df  Sum Sq Mean Sq F value   Pr(>F)   
Age_num     1   797.1  797.06  7.1276 0.008523 **
Urban       1   515.7  515.66  4.6112 0.033550 * 
Residuals 135 15096.7  111.83                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Age_num + Urban"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-17.219  -8.348  -1.427   5.631  30.107 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 14.79556    1.66060   8.910 3.08e-15 ***
Age_num      0.16844    0.05563   3.028  0.00295 ** 
UrbanMedium  4.11374    1.91571   2.147  0.03355 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 10.57 on 135 degrees of freedom
Multiple R-squared:   0.08, Adjusted R-squared:  0.06637 
F-statistic: 5.869 on 2 and 135 DF,  p-value: 0.003595

Next we perform the broken stick regression based on the selected final model

stick_models <- list()
for (YY in Yrs) {
    # YY = '2015'
    print(paste0("-------Following is for ", YY, "--------------"))
    Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
        Age_num >= 1)
    fit.seg <- segmented(lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
    summary(fit.seg) %>% print()
    plot(fit.seg$fitted.values, fit.seg$residuals)
    plot(fit.seg$model$faith_pd, fit.seg$fitted.values)
    qqnorm(fit.seg$residuals)
    if (!"segmented" %in% class(fit.seg)) {
        print(paste(YY, BB, AA, "same as lm model", sep = "-"))
    } else {
        anova(fit.seg, lm_models[[YY]]) %>% print()
        davies.test(lm_models[[YY]], seg.Z = ~Age_num, k = 10) %>% print()
        confint.segmented(fit.seg) %>% print()
        stick_models[[YY]] <- fit.seg
        slope(fit.seg) %>% print()
    }
    
    
}
[1] "-------Following is for 2015--------------"

    ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)

Estimated Break-Point(s):
                Est. St.Err
psi1.Age_num 10.503  1.915

Meaningful coefficients of the linear terms:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   1.2593     3.5611   0.354 0.724183    
Age_num       2.0584     0.5730   3.592 0.000461 ***
UrbanMedium   7.4982     2.1762   3.446 0.000764 ***
U1.Age_num   -2.1586     0.5797  -3.724       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 12.12 on 132 degrees of freedom
Multiple R-Squared: 0.2359,  Adjusted R-squared: 0.2127 

Convergence attained in 3 iter. (rel. change 0)

Analysis of Variance Table

Model 1: faith_pd ~ Age_num + Urban + U1.Age_num + psi1.Age_num
Model 2: faith_pd ~ Age_num + Urban
  Res.Df   RSS Df Sum of Sq      F    Pr(>F)    
1    132 19395                                  
2    134 22985 -2   -3589.7 12.216 1.357e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    Davies' test for a change in the slope

data:  formula = faith_pd ~ Age_num + Urban ,   method = lm 
model = gaussian , link = identity  
segmented variable = Age_num
'best' at = 14.111, n.points = 9, p-value = 7.96e-05
alternative hypothesis: two.sided

                Est. CI(95%).low CI(95%).up
psi1.Age_num 10.5028     6.71477    14.2907
$Age_num
           Est.  St.Err. t value CI(95%).l CI(95%).u
slope1  2.05840 0.573020  3.5922   0.92492  3.191900
slope2 -0.10021 0.096118 -1.0426  -0.29034  0.089921

[1] "-------Following is for 2016--------------"

    ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)

Estimated Break-Point(s):
              Est. St.Err
psi1.Age_num   13  2.234

Meaningful coefficients of the linear terms:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   5.3916     2.7305   1.975   0.0504 .  
Age_num       1.3496     0.3351   4.028 9.41e-05 ***
UrbanMedium   4.6226     1.7888   2.584   0.0108 *  
U1.Age_num   -1.4254     0.3456  -4.125       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.854 on 133 degrees of freedom
Multiple R-Squared: 0.213,  Adjusted R-squared: 0.1894 

Convergence attained in 1 iter. (rel. change -2.1389e-07)

Analysis of Variance Table

Model 1: faith_pd ~ Age_num + Urban + U1.Age_num + psi1.Age_num
Model 2: faith_pd ~ Age_num + Urban
  Res.Df   RSS Df Sum of Sq      F    Pr(>F)    
1    133 12914                                  
2    135 15097 -2   -2183.2 11.242 3.082e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    Davies' test for a change in the slope

data:  formula = faith_pd ~ Age_num + Urban ,   method = lm 
model = gaussian , link = identity  
segmented variable = Age_num
'best' at = 15.444, n.points = 9, p-value = 0.0002563
alternative hypothesis: two.sided

                Est. CI(95%).low CI(95%).up
psi1.Age_num 12.9999     8.58038    17.4194
$Age_num
            Est.  St.Err.  t value CI(95%).l CI(95%).u
slope1  1.349600 0.335070  4.02800   0.68689  2.012400
slope2 -0.075798 0.087213 -0.86912  -0.24830  0.096705

Stick model is better than linear

data for plots

dat15 <- Work %>% filter(Body_Site == BB, Year == "2015", SampleGroup == "Villagers", 
    Age_num >= 1)
dat16 <- Work %>% filter(Body_Site == BB, Year == "2016", SampleGroup == "Villagers", 
    Age_num >= 1)

dat15 <- cbind(dat15, predict(stick_models[["2015"]], newdata = dat15, interval = "confidence"))
dat16 <- cbind(dat16, predict(stick_models[["2016"]], newdata = dat16, interval = "confidence"))

# lm model children between 1-8
for (YY in Yrs) {
    Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
        Age_num >= 1, Age_num <= 8)
    fit.child <- lm(paste0(AA, " ~ Age_num + Urban"), data = Work_dat)
    print(paste0("--------Below is ANOVA table for year ", YY))
    anova(fit.child) %>% print()
    
    print(paste0("--------Below is coefficient for year ", YY))
    summary(fit.child) %>% print()
}
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: faith_pd
          Df Sum Sq Mean Sq F value  Pr(>F)  
Age_num    1  659.2  659.21  3.8258 0.05643 .
Urban      1  300.2  300.19  1.7422 0.19326  
Residuals 47 8098.4  172.31                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Age_num + Urban"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-14.891  -7.674  -3.182   2.997  41.943 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)   3.1387     4.7677   0.658   0.5135  
Age_num       1.8759     0.8572   2.189   0.0336 *
UrbanMedium   5.0475     3.8241   1.320   0.1933  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 13.13 on 47 degrees of freedom
Multiple R-squared:  0.1059,    Adjusted R-squared:  0.06787 
F-statistic: 2.784 on 2 and 47 DF,  p-value: 0.072

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: faith_pd
          Df Sum Sq Mean Sq F value  Pr(>F)  
Age_num    1  121.3  121.33  1.1391 0.29166  
Urban      1  415.6  415.63  3.9020 0.05452 .
Residuals 44 4686.8  106.52                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Age_num + Urban"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-13.915  -6.450  -2.914   3.508  31.106 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)   6.4663     4.1339   1.564   0.1249  
Age_num       0.9505     0.7221   1.316   0.1949  
UrbanMedium   6.0112     3.0431   1.975   0.0545 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 10.32 on 44 degrees of freedom
Multiple R-squared:  0.1028,    Adjusted R-squared:  0.06201 
F-statistic: 2.521 on 2 and 44 DF,  p-value: 0.09197
dat_p[["2015"]][[AA]] = dat15
dat_p[["2016"]][[AA]] = dat16
observed_otu
BB = "Nose"
AA = "observed_otus"

Work_dat <- Work %>% filter(Body_Site == BB, SampleGroup == "Villagers", Age_num >= 
    1)
p <- ggplot(Work_dat, aes(x = Age_num, y = get(AA), color = Urban, shape = Gender)) + 
    geom_point(size = 1.5) + scale_shape_manual(values = c(16, 21)) + scale_color_manual(breaks = names(Urban_color), 
    values = Urban_color) + facet_grid(. ~ Year) + theme_jw() + theme(aspect.ratio = 0.8, 
    panel.background = element_rect(color = "black")) + labs(x = "Age (Years)", y = AA, 
    title = BB)
print(p)

# Year 2015
YY = "2015"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
    Age_num >= 1)

# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender + 
    Urban:Gender, data = Work_dat)

# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)

# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num + 
    Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)

anova(regforward.out)
# Year 2016
YY = "2016"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
    Age_num >= 1)

# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender + 
    Urban:Gender, data = Work_dat)

# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)

# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num + 
    Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)

anova(regforward.out)
# Final linear model
lm_models <- list()
for (YY in Yrs) {
    Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
        Age_num >= 1)
    fit.final <- lm(paste0(AA, " ~ Age_num + Urban"), data = Work_dat)
    print(paste0("--------Below is ANOVA table for year ", YY))
    anova(fit.final) %>% print()
    
    print(paste0("--------Below is coefficient for year ", YY))
    summary(fit.final) %>% print()
    lm_models[[YY]] <- fit.final
}
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: observed_otus
           Df  Sum Sq Mean Sq F value  Pr(>F)  
Age_num     1   19166   19166  1.9738 0.16236  
Urban       1   57605   57605  5.9325 0.01618 *
Residuals 134 1301147    9710                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Age_num + Urban"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-113.75  -61.06  -27.67   38.30  433.79 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  76.9081    15.3569   5.008  1.7e-06 ***
Age_num       0.9394     0.5202   1.806   0.0732 .  
UrbanMedium  42.9052    17.6153   2.436   0.0162 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 98.54 on 134 degrees of freedom
Multiple R-squared:  0.05571,   Adjusted R-squared:  0.04162 
F-statistic: 3.953 on 2 and 134 DF,  p-value: 0.02147

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: observed_otus
           Df Sum Sq Mean Sq F value  Pr(>F)  
Age_num     1  24477 24477.4  4.6338 0.03313 *
Urban       1  10229 10229.2  1.9365 0.16634  
Residuals 135 713122  5282.4                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Age_num + Urban"), data = Work_dat)

Residuals:
   Min     1Q Median     3Q    Max 
-90.69 -54.83 -16.74  27.42 270.39 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  76.7036    11.4132   6.721 4.64e-10 ***
Age_num       0.9089     0.3824   2.377   0.0189 *  
UrbanMedium  18.3221    13.1665   1.392   0.1663    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 72.68 on 135 degrees of freedom
Multiple R-squared:  0.04641,   Adjusted R-squared:  0.03228 
F-statistic: 3.285 on 2 and 135 DF,  p-value: 0.04045

Next we perform the broken stick regression based on the selected final model

stick_models <- list()
for (YY in Yrs) {
    # YY = '2015'
    Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
        Age_num >= 1)
    fit.seg <- segmented(lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
    summary(fit.seg) %>% print()
    plot(fit.seg$fitted.values, fit.seg$residuals)
    plot(fit.seg$model$observed_otus, fit.seg$fitted.values)
    qqnorm(fit.seg$residuals)
    
    if (!"segmented" %in% class(fit.seg)) {
        print(paste(YY, BB, AA, "same as lm model", sep = "-"))
    } else {
        anova(fit.seg, lm_models[[YY]]) %>% print()
        davies.test(lm_models[[YY]], seg.Z = ~Age_num, k = 10) %>% print()
        confint.segmented(fit.seg) %>% print()
        stick_models[[YY]] <- fit.seg
        slope(fit.seg) %>% print()
    }
    
    
}

    ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)

Estimated Break-Point(s):
              Est. St.Err
psi1.Age_num   10   1.86

Meaningful coefficients of the linear terms:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -14.592     27.044  -0.540 0.590405    
Age_num       15.130      4.352   3.477 0.000687 ***
UrbanMedium   45.538     16.526   2.755 0.006690 ** 
U1.Age_num   -16.118      4.402  -3.661       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 92.05 on 132 degrees of freedom
Multiple R-Squared: 0.1882,  Adjusted R-squared: 0.1636 

Convergence attained in 6 iter. (rel. change 1.8623e-07)

Analysis of Variance Table

Model 1: observed_otus ~ Age_num + Urban + U1.Age_num + psi1.Age_num
Model 2: observed_otus ~ Age_num + Urban
  Res.Df     RSS Df Sum of Sq      F    Pr(>F)    
1    132 1118559                                  
2    134 1301147 -2   -182588 10.774 4.634e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    Davies' test for a change in the slope

data:  formula = observed_otus ~ Age_num + Urban ,   method = lm 
model = gaussian , link = identity  
segmented variable = Age_num
'best' at = 7.5556, n.points = 9, p-value = 0.0002517
alternative hypothesis: two.sided

             Est. CI(95%).low CI(95%).up
psi1.Age_num   10      6.3204    13.6796
$Age_num
         Est. St.Err. t value CI(95%).l CI(95%).u
slope1 15.130 4.35170  3.4769    6.5225   23.7380
slope2 -0.988 0.72994 -1.3535   -2.4319    0.4559


    ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)

Estimated Break-Point(s):
              Est. St.Err
psi1.Age_num   13  2.435

Meaningful coefficients of the linear terms:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   16.708     18.984   0.880 0.380400    
Age_num        8.445      2.330   3.625 0.000411 ***
UrbanMedium   21.568     12.437   1.734 0.085200 .  
U1.Age_num    -9.094      2.403  -3.785       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 68.51 on 133 degrees of freedom
Multiple R-Squared: 0.1652,  Adjusted R-squared: 0.1401 

Convergence attained in 4 iter. (rel. change 1.796e-07)

Analysis of Variance Table

Model 1: observed_otus ~ Age_num + Urban + U1.Age_num + psi1.Age_num
Model 2: observed_otus ~ Age_num + Urban
  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1    133 624263                                  
2    135 713122 -2    -88860 9.4658 0.0001434 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    Davies' test for a change in the slope

data:  formula = observed_otus ~ Age_num + Urban ,   method = lm 
model = gaussian , link = identity  
segmented variable = Age_num
'best' at = 15.444, n.points = 9, p-value = 0.000653
alternative hypothesis: two.sided

             Est. CI(95%).low CI(95%).up
psi1.Age_num   13     8.18352    17.8164
$Age_num
           Est. St.Err. t value CI(95%).l CI(95%).u
slope1  8.44460 2.32960  3.6249    3.8367  13.05300
slope2 -0.64936 0.60637 -1.0709   -1.8487   0.55002

data for plots

dat15 <- Work %>% filter(Body_Site == BB, Year == "2015", SampleGroup == "Villagers", 
    Age_num >= 1)
dat16 <- Work %>% filter(Body_Site == BB, Year == "2016", SampleGroup == "Villagers", 
    Age_num >= 1)

dat15 <- cbind(dat15, predict(stick_models[["2015"]], newdata = dat15, interval = "confidence"))
dat16 <- cbind(dat16, predict(stick_models[["2016"]], newdata = dat16, interval = "confidence"))

# lm model children
for (YY in Yrs) {
    Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
        Age_num >= 1, Age_num <= 8)
    fit.child <- lm(paste0(AA, " ~ Age_num + Urban"), data = Work_dat)
    print(paste0("--------Below is ANOVA table for year ", YY))
    anova(fit.child) %>% print()
    
    print(paste0("--------Below is coefficient for year ", YY))
    summary(fit.child) %>% print()
}
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: observed_otus
          Df Sum Sq Mean Sq F value  Pr(>F)  
Age_num    1  35060   35060  3.9447 0.05287 .
Urban      1  10828   10828  1.2183 0.27531  
Residuals 47 417726    8888                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Age_num + Urban"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-115.57  -54.57  -29.85   13.78  298.43 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  -0.7517    34.2415  -0.022   0.9826  
Age_num      13.3763     6.1561   2.173   0.0349 *
UrbanMedium  30.3149    27.4648   1.104   0.2753  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 94.28 on 47 degrees of freedom
Multiple R-squared:  0.09898,   Adjusted R-squared:  0.06064 
F-statistic: 2.582 on 2 and 47 DF,  p-value: 0.08635

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: observed_otus
          Df Sum Sq Mean Sq F value Pr(>F)
Age_num    1    737   737.0  0.1750 0.6777
Urban      1   5715  5715.4  1.3571 0.2503
Residuals 44 185304  4211.5               
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Age_num + Urban"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-66.219 -36.428 -23.078   6.143 207.781 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   41.630     25.993   1.602    0.116
Age_num        2.575      4.541   0.567    0.574
UrbanMedium   22.291     19.135   1.165    0.250

Residual standard error: 64.9 on 44 degrees of freedom
Multiple R-squared:  0.03365,   Adjusted R-squared:  -0.01028 
F-statistic: 0.7661 on 2 and 44 DF,  p-value: 0.4709
dat_p[["2015"]][[AA]] = dat15
dat_p[["2016"]][[AA]] = dat16
shannon
BB = "Nose"
AA = "shannon"

Work_dat <- Work %>% filter(Body_Site == BB, SampleGroup == "Villagers", Age_num >= 
    1)
p <- ggplot(Work_dat, aes(x = Age_num, y = get(AA), color = Urban, shape = Gender)) + 
    geom_point(size = 1.5) + scale_shape_manual(values = c(16, 21)) + scale_color_manual(breaks = names(Urban_color), 
    values = Urban_color) + facet_grid(. ~ Year) + theme_jw() + theme(aspect.ratio = 0.8, 
    panel.background = element_rect(color = "black")) + labs(x = "Age (Years)", y = AA, 
    title = BB)
print(p)

# Year 2015
YY = "2015"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
    Age_num >= 1)

# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender + 
    Urban:Gender, data = Work_dat)

# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)

# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num + 
    Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)

anova(regforward.out)
# Year 2016
YY = "2016"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
    Age_num >= 1)

# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender + 
    Urban:Gender, data = Work_dat)

# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)

# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num + 
    Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)

anova(regforward.out)

Urban is not significant in both years

# Final linear model
lm_models <- list()
for (YY in Yrs) {
    Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
        Age_num >= 1)
    fit.final <- lm(paste0(AA, " ~ Age_num"), data = Work_dat)
    print(paste0("--------Below is ANOVA table for year ", YY))
    anova(fit.final) %>% print()
    
    print(paste0("--------Below is coefficient for year ", YY))
    summary(fit.final) %>% print()
    lm_models[[YY]] <- fit.final
}
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: shannon
           Df  Sum Sq Mean Sq F value   Pr(>F)   
Age_num     1  11.523 11.5234  7.4168 0.007316 **
Residuals 135 209.750  1.5537                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Age_num"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.4027 -0.7720 -0.1378  0.6977  3.9535 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.84986    0.16264  17.522  < 2e-16 ***
Age_num      0.01765    0.00648   2.723  0.00732 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.246 on 135 degrees of freedom
Multiple R-squared:  0.05208,   Adjusted R-squared:  0.04506 
F-statistic: 7.417 on 1 and 135 DF,  p-value: 0.007316

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: shannon
           Df  Sum Sq Mean Sq F value  Pr(>F)  
Age_num     1   7.564  7.5644  5.8903 0.01653 *
Residuals 136 174.653  1.2842                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Age_num"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.4861 -0.7870 -0.1429  0.5431  3.9504 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 2.893257   0.150401  19.237   <2e-16 ***
Age_num     0.014208   0.005854   2.427   0.0165 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.133 on 136 degrees of freedom
Multiple R-squared:  0.04151,   Adjusted R-squared:  0.03447 
F-statistic:  5.89 on 1 and 136 DF,  p-value: 0.01653

Next we perform the broken stick regression based on the selected final model

stick_models <- list()
for (YY in Yrs) {
    # YY = '2015'
    Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
        Age_num >= 1)
    fit.seg <- segmented(lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
    summary(fit.seg) %>% print()
    plot(fit.seg$fitted.values, fit.seg$residuals)
    plot(fit.seg$model$shannon, fit.seg$fitted.values)
    qqnorm(fit.seg$residuals)
    
    if (!"segmented" %in% class(fit.seg)) {
        print(paste(YY, BB, AA, "same as lm model", sep = "-"))
    } else {
        anova(fit.seg, lm_models[[YY]]) %>% print()
        davies.test(lm_models[[YY]], seg.Z = ~Age_num, k = 10) %>% print()
        confint.segmented(fit.seg) %>% print()
        stick_models[[YY]] <- fit.seg
        slope(fit.seg) %>% print()
    }
    
    
}

    ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)

Estimated Break-Point(s):
              Est. St.Err
psi1.Age_num   12   2.51

Meaningful coefficients of the linear terms:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.94561    0.29349   6.629  7.7e-10 ***
Age_num      0.14682    0.04321   3.398 0.000897 ***
U1.Age_num  -0.15485    0.04441  -3.487       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.182 on 133 degrees of freedom
Multiple R-Squared: 0.1598,  Adjusted R-squared: 0.1409 

Convergence attained in 1 iter. (rel. change 9.1184e-09)

Analysis of Variance Table

Model 1: shannon ~ Age_num + U1.Age_num + psi1.Age_num
Model 2: shannon ~ Age_num
  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1    133 185.90                                  
2    135 209.75 -2   -23.847 8.5303 0.0003269 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    Davies' test for a change in the slope

data:  formula = shannon ~ Age_num ,   method = lm 
model = gaussian , link = identity  
segmented variable = Age_num
'best' at = 14.111, n.points = 9, p-value = 0.0006814
alternative hypothesis: two.sided

                Est. CI(95%).low CI(95%).up
psi1.Age_num 12.0005     7.03577    16.9652
$Age_num
             Est.  St.Err.  t value CI(95%).l CI(95%).u
slope1  0.1468200 0.043210  3.39770  0.061349   0.23228
slope2 -0.0080328 0.010259 -0.78297 -0.028325   0.01226


    ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)

Estimated Break-Point(s):
                Est. St.Err
psi1.Age_num 14.965  2.962

Meaningful coefficients of the linear terms:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.06595    0.27304   7.567  5.5e-12 ***
Age_num      0.11406    0.03512   3.248  0.00147 ** 
U1.Age_num  -0.12816    0.03639  -3.521       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.067 on 134 degrees of freedom
Multiple R-Squared: 0.1627,  Adjusted R-squared: 0.1439 

Convergence attained in 3 iter. (rel. change 2.7832e-07)

Analysis of Variance Table

Model 1: shannon ~ Age_num + U1.Age_num + psi1.Age_num
Model 2: shannon ~ Age_num
  Res.Df    RSS Df Sum of Sq      F   Pr(>F)    
1    134 152.58                                 
2    136 174.65 -2   -22.074 9.6931 0.000117 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    Davies' test for a change in the slope

data:  formula = shannon ~ Age_num ,   method = lm 
model = gaussian , link = identity  
segmented variable = Age_num
'best' at = 15.444, n.points = 9, p-value = 0.0002751
alternative hypothesis: two.sided

               Est. CI(95%).low CI(95%).up
psi1.Age_num 14.965      9.1063    20.8238
$Age_num
            Est.   St.Err. t value CI(95%).l CI(95%).u
slope1  0.114060 0.0351210  3.2476  0.044596 0.1835200
slope2 -0.014103 0.0095434 -1.4777 -0.032978 0.0047725

data for plots

dat15 <- Work %>% filter(Body_Site == BB, Year == "2015", SampleGroup == "Villagers", 
    Age_num >= 1)
dat16 <- Work %>% filter(Body_Site == BB, Year == "2016", SampleGroup == "Villagers", 
    Age_num >= 1)

dat15 <- cbind(dat15, predict(stick_models[["2015"]], newdata = dat15, interval = "confidence"))
dat16 <- cbind(dat16, predict(stick_models[["2016"]], newdata = dat16, interval = "confidence"))

dat_p[["2015"]][[AA]] = dat15
dat_p[["2016"]][[AA]] = dat16
# lm model children
for (YY in Yrs) {
    Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
        Age_num >= 1, Age_num <= 8)
    fit.child <- lm(paste0(AA, " ~ Age_num + Urban"), data = Work_dat)
    print(paste0("--------Below is ANOVA table for year ", YY))
    anova(fit.child) %>% print()
    
    print(paste0("--------Below is coefficient for year ", YY))
    summary(fit.child) %>% print()
}
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: shannon
          Df Sum Sq Mean Sq F value Pr(>F)
Age_num    1  3.638  3.6384  2.7025 0.1069
Urban      1  1.025  1.0249  0.7613 0.3874
Residuals 47 63.275  1.3463               
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Age_num + Urban"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.4719 -0.6517 -0.2165  0.4322  4.5073 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.22962    0.42143   5.291 3.13e-06 ***
Age_num      0.10800    0.07577   1.425    0.161    
UrbanMedium -0.29494    0.33802  -0.873    0.387    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.16 on 47 degrees of freedom
Multiple R-squared:  0.06864,   Adjusted R-squared:  0.02901 
F-statistic: 1.732 on 2 and 47 DF,  p-value: 0.188

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: shannon
          Df Sum Sq Mean Sq F value Pr(>F)
Age_num    1  0.228 0.22845  0.2025 0.6550
Urban      1  0.979 0.97901  0.8676 0.3567
Residuals 44 49.649 1.12840               
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Age_num + Urban"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.9722 -0.5915 -0.1220  0.2949  2.9743 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.89677    0.42548   6.808 2.18e-08 ***
Age_num     -0.04220    0.07432  -0.568    0.573    
UrbanMedium -0.29175    0.31321  -0.931    0.357    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.062 on 44 degrees of freedom
Multiple R-squared:  0.02374,   Adjusted R-squared:  -0.02063 
F-statistic: 0.535 on 2 and 44 DF,  p-value: 0.5894
Data for nose
dat_all[[BB]] = dat_p

Right_Arm

dat_p = list()
dat_p[["2015"]] = list()
dat_p[["2016"]] = list()

Modeling without babies

faith_pd
BB = "Right_Arm"
AA = "faith_pd"

Work_dat <- Work %>% filter(Body_Site == BB, SampleGroup == "Villagers", Age_num >= 
    1)
p <- ggplot(Work_dat, aes(x = Age_num, y = get(AA), color = Urban, shape = Gender)) + 
    geom_point(size = 1.5) + scale_shape_manual(values = c(16, 21)) + scale_color_manual(breaks = names(Urban_color), 
    values = Urban_color) + facet_grid(. ~ Year) + theme_jw() + theme(aspect.ratio = 0.8, 
    panel.background = element_rect(color = "black")) + labs(x = "Age (Years)", y = AA, 
    title = BB)
print(p)

# Year 2015
YY = "2015"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
    Age_num >= 1)

# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender + 
    Urban:Gender, data = Work_dat)

# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)

# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num + 
    Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)

anova(regforward.out)

For year 2015, only urban and gender is significant in determing the faith_pd

# Year 2016
YY = "2016"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
    Age_num >= 1)

# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender + 
    Urban:Gender, data = Work_dat)

# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)

# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num + 
    Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)

anova(regforward.out)

In 2016, Only Urban

For consistency, only gender and urban were included.

# Final linear model
lm_models <- list()
for (YY in Yrs) {
    Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
        Age_num >= 1)
    fit.final <- lm(paste0(AA, " ~ Urban + Gender"), data = Work_dat)
    print(paste0("--------Below is ANOVA table for year ", YY))
    anova(fit.final) %>% print()
    
    print(paste0("--------Below is coefficient for year ", YY))
    summary(fit.final) %>% print()
    lm_models[[YY]] <- fit.final
}
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: faith_pd
           Df  Sum Sq Mean Sq F value    Pr(>F)    
Urban       1  9453.0  9453.0  53.661 1.083e-11 ***
Gender      1   978.2   978.2   5.553   0.01965 *  
Residuals 161 28361.9   176.2                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Urban + Gender"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-31.938  -8.529  -0.744   7.129  43.738 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   41.131      1.681  24.470  < 2e-16 ***
UrbanMedium   15.519      2.093   7.413 6.62e-12 ***
GenderM        4.918      2.087   2.356   0.0197 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 13.27 on 161 degrees of freedom
Multiple R-squared:  0.2689,    Adjusted R-squared:  0.2598 
F-statistic: 29.61 on 2 and 161 DF,  p-value: 1.123e-11

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: faith_pd
          Df  Sum Sq Mean Sq F value  Pr(>F)  
Urban      1  1468.4 1468.41  6.3177 0.01396 *
Gender     1    52.2   52.23  0.2247 0.63676  
Residuals 80 18594.1  232.43                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Urban + Gender"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-29.123 -11.062  -2.842  11.114  42.245 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   42.860      2.697  15.895   <2e-16 ***
UrbanMedium    8.388      3.347   2.506   0.0142 *  
GenderM        1.640      3.460   0.474   0.6368    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 15.25 on 80 degrees of freedom
Multiple R-squared:  0.0756,    Adjusted R-squared:  0.05249 
F-statistic: 3.271 on 2 and 80 DF,  p-value: 0.0431

data for plots

dat15 <- Work %>% filter(Body_Site == BB, Year == "2015", SampleGroup == "Villagers", 
    Age_num >= 1)
dat16 <- Work %>% filter(Body_Site == BB, Year == "2016", SampleGroup == "Villagers", 
    Age_num >= 1)

dat15 <- cbind(dat15, predict(lm_models[["2015"]], newdata = dat15, interval = "confidence"))
dat16 <- cbind(dat16, predict(lm_models[["2016"]], newdata = dat16, interval = "confidence"))

dat_p[["2015"]][[AA]] = dat15
dat_p[["2016"]][[AA]] = dat16
observed_otu
BB = "Right_Arm"
AA = "observed_otus"

Work_dat <- Work %>% filter(Body_Site == BB, SampleGroup == "Villagers", Age_num >= 
    1)
p <- ggplot(Work_dat, aes(x = Age_num, y = get(AA), color = Urban, shape = Gender)) + 
    geom_point(size = 1.5) + scale_shape_manual(values = c(16, 21)) + scale_color_manual(breaks = names(Urban_color), 
    values = Urban_color) + facet_grid(. ~ Year) + theme_jw() + theme(aspect.ratio = 0.8, 
    panel.background = element_rect(color = "black")) + labs(x = "Age (Years)", y = AA, 
    title = BB)
print(p)

# Year 2015
YY = "2015"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
    Age_num >= 1)

# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender + 
    Urban:Gender, data = Work_dat)

# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)

# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num + 
    Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)

anova(regforward.out)
# Year 2016
YY = "2016"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
    Age_num >= 1)

# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender + 
    Urban:Gender, data = Work_dat)

# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)

# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num + 
    Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)

anova(regforward.out)
# Final linear model
lm_models <- list()
for (YY in Yrs) {
    Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
        Age_num >= 1)
    fit.final <- lm(paste0(AA, " ~ Urban + Gender"), data = Work_dat)
    print(paste0("--------Below is ANOVA table for year ", YY))
    anova(fit.final) %>% print()
    
    print(paste0("--------Below is coefficient for year ", YY))
    summary(fit.final) %>% print()
    lm_models[[YY]] <- fit.final
}
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: observed_otus
           Df  Sum Sq Mean Sq F value    Pr(>F)    
Urban       1  415511  415511 34.9526 1.959e-08 ***
Gender      1   71908   71908  6.0488   0.01497 *  
Residuals 161 1913942   11888                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Urban + Gender"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-243.34  -71.27  -10.54   53.40  350.66 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   300.90      13.81  21.791  < 2e-16 ***
UrbanMedium   103.27      17.20   6.005 1.23e-08 ***
GenderM        42.17      17.14   2.459    0.015 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 109 on 161 degrees of freedom
Multiple R-squared:  0.203, Adjusted R-squared:  0.1931 
F-statistic:  20.5 on 2 and 161 DF,  p-value: 1.171e-08

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: observed_otus
          Df  Sum Sq Mean Sq F value  Pr(>F)  
Urban      1   62837   62837  3.2078 0.07707 .
Gender     1    4018    4018  0.2051 0.65185  
Residuals 80 1567098   19589                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Urban + Gender"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-235.71 -114.92  -30.32  101.79  454.29 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   326.32      24.75  13.182   <2e-16 ***
UrbanMedium    54.82      30.73   1.784   0.0783 .  
GenderM        14.39      31.76   0.453   0.6518    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 140 on 80 degrees of freedom
Multiple R-squared:  0.04092,   Adjusted R-squared:  0.01694 
F-statistic: 1.706 on 2 and 80 DF,  p-value: 0.188

data for plots

dat15 <- Work %>% filter(Body_Site == BB, Year == "2015", SampleGroup == "Villagers", 
    Age_num >= 1)
dat16 <- Work %>% filter(Body_Site == BB, Year == "2016", SampleGroup == "Villagers", 
    Age_num >= 1)

dat15 <- cbind(dat15, predict(lm_models[["2015"]], newdata = dat15, interval = "confidence"))
dat16 <- cbind(dat16, predict(lm_models[["2016"]], newdata = dat16, interval = "confidence"))

dat_p[["2015"]][[AA]] = dat15
dat_p[["2016"]][[AA]] = dat16
shannon
BB = "Right_Arm"
AA = "shannon"

Work_dat <- Work %>% filter(Body_Site == BB, SampleGroup == "Villagers", Age_num >= 
    1)
p <- ggplot(Work_dat, aes(x = Age_num, y = get(AA), color = Urban, shape = Gender)) + 
    geom_point(size = 1.5) + scale_shape_manual(values = c(16, 21)) + scale_color_manual(breaks = names(Urban_color), 
    values = Urban_color) + facet_grid(. ~ Year) + theme_jw() + theme(aspect.ratio = 0.8, 
    panel.background = element_rect(color = "black")) + labs(x = "Age (Years)", y = AA, 
    title = BB)
print(p)

# Year 2015
YY = "2015"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
    Age_num >= 1)

# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender + 
    Urban:Gender, data = Work_dat)

# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)

# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num + 
    Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)

anova(regforward.out)
# Year 2016
YY = "2016"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
    Age_num >= 1)

# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender + 
    Urban:Gender, data = Work_dat)

# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)

# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num + 
    Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)

anova(regforward.out)

Urban is not significant in both years

# Final linear model
lm_models <- list()
for (YY in Yrs) {
    Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
        Age_num >= 1)
    fit.final <- lm(paste0(AA, " ~ Urban + Gender"), data = Work_dat)
    print(paste0("--------Below is ANOVA table for year ", YY))
    anova(fit.final) %>% print()
    
    print(paste0("--------Below is coefficient for year ", YY))
    summary(fit.final) %>% print()
    lm_models[[YY]] <- fit.final
}
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: shannon
           Df Sum Sq Mean Sq F value    Pr(>F)    
Urban       1  6.697  6.6969 11.3078 0.0009644 ***
Gender      1  5.370  5.3699  9.0671 0.0030226 ** 
Residuals 161 95.350  0.5922                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Urban + Gender"), data = Work_dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.90249 -0.44609  0.04411  0.51664  1.58334 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  6.36287    0.09746  65.286  < 2e-16 ***
UrbanMedium  0.42235    0.12138   3.480 0.000646 ***
GenderM      0.36439    0.12101   3.011 0.003023 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7696 on 161 degrees of freedom
Multiple R-squared:  0.1123,    Adjusted R-squared:  0.1013 
F-statistic: 10.19 on 2 and 161 DF,  p-value: 6.824e-05

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: shannon
          Df Sum Sq Mean Sq F value   Pr(>F)    
Urban      1  6.884  6.8840 11.6908 0.000991 ***
Gender     1  0.028  0.0284  0.0482 0.826864    
Residuals 80 47.108  0.5888                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Urban + Gender"), data = Work_dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.99333 -0.34828 -0.00307  0.56205  1.52427 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  6.44814    0.13572  47.509  < 2e-16 ***
UrbanMedium  0.57661    0.16849   3.422 0.000981 ***
GenderM     -0.03821    0.17414  -0.219 0.826864    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7674 on 80 degrees of freedom
Multiple R-squared:  0.128, Adjusted R-squared:  0.1062 
F-statistic: 5.869 on 2 and 80 DF,  p-value: 0.004183

data for plots

dat15 <- Work %>% filter(Body_Site == BB, Year == "2015", SampleGroup == "Villagers", 
    Age_num >= 1)
dat16 <- Work %>% filter(Body_Site == BB, Year == "2016", SampleGroup == "Villagers", 
    Age_num >= 1)

dat15 <- cbind(dat15, predict(lm_models[["2015"]], newdata = dat15, interval = "confidence"))
dat16 <- cbind(dat16, predict(lm_models[["2016"]], newdata = dat16, interval = "confidence"))

dat_p[["2015"]][[AA]] = dat15
dat_p[["2016"]][[AA]] = dat16
Data for right arm
dat_all[[BB]] = dat_p

Right_Hand

dat_p = list()
dat_p[["2015"]] = list()
dat_p[["2016"]] = list()

Modeling without babies

faith_pd
BB = "Right_Hand"
AA = "faith_pd"

Work_dat <- Work %>% filter(Body_Site == BB, SampleGroup == "Villagers", Age_num >= 
    1)
p <- ggplot(Work_dat, aes(x = Age_num, y = get(AA), color = Urban, shape = Gender)) + 
    geom_point(size = 1.5) + scale_shape_manual(values = c(16, 21)) + scale_color_manual(breaks = names(Urban_color), 
    values = Urban_color) + facet_grid(. ~ Year) + theme_jw() + theme(aspect.ratio = 0.8, 
    panel.background = element_rect(color = "black")) + labs(x = "Age (Years)", y = AA, 
    title = BB)
print(p)

# Year 2015
YY = "2015"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
    Age_num >= 1)

# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender + 
    Urban:Gender, data = Work_dat)

# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)

# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num + 
    Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)

anova(regforward.out)
# Year 2016
YY = "2016"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
    Age_num >= 1)

# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender + 
    Urban:Gender, data = Work_dat)

# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)

# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num + 
    Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)

anova(regforward.out)

In 2016, Only Urban

For consistency, only gender and urban were included.

# Final linear model
lm_models <- list()
for (YY in Yrs) {
    Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
        Age_num >= 1)
    fit.final <- lm(paste0(AA, " ~ Urban + Gender + Age_num"), data = Work_dat)
    print(paste0("--------Below is ANOVA table for year ", YY))
    anova(fit.final) %>% print()
    
    print(paste0("--------Below is coefficient for year ", YY))
    summary(fit.final) %>% print()
    lm_models[[YY]] <- fit.final
}
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: faith_pd
           Df  Sum Sq Mean Sq F value    Pr(>F)    
Urban       1 17233.1 17233.1 121.850 < 2.2e-16 ***
Gender      1  1705.5  1705.5  12.059 0.0006631 ***
Age_num     1   588.1   588.1   4.158 0.0430834 *  
Residuals 160 22628.6   141.4                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Urban + Gender + Age_num"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-35.668  -7.504   0.730   6.352  33.437 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 39.10624    1.91115  20.462  < 2e-16 ***
UrbanMedium 20.24935    1.89466  10.688  < 2e-16 ***
GenderM      6.61860    1.87575   3.529 0.000546 ***
Age_num     -0.11885    0.05829  -2.039 0.043083 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 11.89 on 160 degrees of freedom
Multiple R-squared:  0.4632,    Adjusted R-squared:  0.4531 
F-statistic: 46.02 on 3 and 160 DF,  p-value: < 2.2e-16

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: faith_pd
           Df  Sum Sq Mean Sq F value    Pr(>F)    
Urban       1  7975.2  7975.2 36.5302 1.978e-08 ***
Gender      1   799.9   799.9  3.6638   0.05813 .  
Age_num     1   545.3   545.3  2.4976   0.11681    
Residuals 113 24669.9   218.3                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Urban + Gender + Age_num"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-34.010  -9.883  -0.960  11.680  29.544 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 37.56772    2.70294  13.899  < 2e-16 ***
UrbanMedium 16.51829    2.74989   6.007 2.35e-08 ***
GenderM      5.43211    2.82161   1.925   0.0567 .  
Age_num     -0.13625    0.08621  -1.580   0.1168    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 14.78 on 113 degrees of freedom
Multiple R-squared:  0.2742,    Adjusted R-squared:  0.2549 
F-statistic: 14.23 on 3 and 113 DF,  p-value: 6.252e-08

Next we perform the broken stick regression based on the selected final model

stick_models <- list()
for (YY in Yrs) {
    # YY = '2015'
    print(paste0("-------Following is for ", YY, "--------------"))
    Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
        Age_num >= 1)
    fit.seg <- segmented(lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
    summary(fit.seg) %>% print()
    plot(fit.seg$fitted.values, fit.seg$residuals)
    plot(fit.seg$model$faith_pd, fit.seg$fitted.values)
    qqnorm(fit.seg$residuals)
    if (!"segmented" %in% class(fit.seg)) {
        print(paste(YY, BB, AA, "same as lm model", sep = "-"))
    } else {
        anova(fit.seg, lm_models[[YY]]) %>% print()
        davies.test(lm_models[[YY]], seg.Z = ~Age_num, k = 10) %>% print()
        confint.segmented(fit.seg) %>% print()
        stick_models[[YY]] <- fit.seg
        slope(fit.seg) %>% print()
    }
}
[1] "-------Following is for 2015--------------"

Call:
lm(formula = paste0(AA, " ~ Urban + Gender + Age_num"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-35.668  -7.504   0.730   6.352  33.437 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 39.10624    1.91115  20.462  < 2e-16 ***
UrbanMedium 20.24935    1.89466  10.688  < 2e-16 ***
GenderM      6.61860    1.87575   3.529 0.000546 ***
Age_num     -0.11885    0.05829  -2.039 0.043083 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 11.89 on 160 degrees of freedom
Multiple R-squared:  0.4632,    Adjusted R-squared:  0.4531 
F-statistic: 46.02 on 3 and 160 DF,  p-value: < 2.2e-16

[1] "2015-Right_Hand-faith_pd-same as lm model"
[1] "-------Following is for 2016--------------"

Call:
lm(formula = paste0(AA, " ~ Urban + Gender + Age_num"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-34.010  -9.883  -0.960  11.680  29.544 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 37.56772    2.70294  13.899  < 2e-16 ***
UrbanMedium 16.51829    2.74989   6.007 2.35e-08 ***
GenderM      5.43211    2.82161   1.925   0.0567 .  
Age_num     -0.13625    0.08621  -1.580   0.1168    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 14.78 on 113 degrees of freedom
Multiple R-squared:  0.2742,    Adjusted R-squared:  0.2549 
F-statistic: 14.23 on 3 and 113 DF,  p-value: 6.252e-08

[1] "2016-Right_Hand-faith_pd-same as lm model"

data for plots

dat15 <- Work %>% filter(Body_Site == BB, Year == "2015", SampleGroup == "Villagers", 
    Age_num >= 1)
dat16 <- Work %>% filter(Body_Site == BB, Year == "2016", SampleGroup == "Villagers", 
    Age_num >= 1)

dat15 <- cbind(dat15, predict(lm_models[["2015"]], newdata = dat15, interval = "confidence"))
dat16 <- cbind(dat16, predict(lm_models[["2016"]], newdata = dat16, interval = "confidence"))

dat_p[["2015"]][[AA]] = dat15
dat_p[["2016"]][[AA]] = dat16
observed_otu
BB = "Right_Hand"
AA = "observed_otus"

Work_dat <- Work %>% filter(Body_Site == BB, SampleGroup == "Villagers", Age_num >= 
    1)
p <- ggplot(Work_dat, aes(x = Age_num, y = get(AA), color = Urban, shape = Gender)) + 
    geom_point(size = 1.5) + scale_shape_manual(values = c(16, 21)) + scale_color_manual(breaks = names(Urban_color), 
    values = Urban_color) + facet_grid(. ~ Year) + theme_jw() + theme(aspect.ratio = 0.8, 
    panel.background = element_rect(color = "black")) + labs(x = "Age (Years)", y = AA, 
    title = BB)
print(p)

# Year 2015
YY = "2015"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
    Age_num >= 1)

# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender + 
    Urban:Gender, data = Work_dat)

# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)

# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num + 
    Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)

anova(regforward.out)
# Year 2016
YY = "2016"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
    Age_num >= 1)

# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender + 
    Urban:Gender, data = Work_dat)

# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)

# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num + 
    Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)

anova(regforward.out)
# Final linear model
lm_models <- list()
for (YY in Yrs) {
    Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
        Age_num >= 1)
    fit.final <- lm(paste0(AA, " ~ Urban + Gender + Age_num"), data = Work_dat)
    print(paste0("--------Below is ANOVA table for year ", YY))
    anova(fit.final) %>% print()
    
    print(paste0("--------Below is coefficient for year ", YY))
    summary(fit.final) %>% print()
    lm_models[[YY]] <- fit.final
}
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: observed_otus
           Df  Sum Sq Mean Sq F value  Pr(>F)    
Urban       1  938567  938567 92.5543 < 2e-16 ***
Gender      1   76969   76969  7.5901 0.00655 ** 
Age_num     1   34657   34657  3.4176 0.06635 .  
Residuals 160 1622514   10141                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Urban + Gender + Age_num"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-253.91  -66.40   -0.44   55.81  354.88 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 286.3689    16.1831  17.696  < 2e-16 ***
UrbanMedium 149.1367    16.0434   9.296  < 2e-16 ***
GenderM      44.5666    15.8833   2.806  0.00564 ** 
Age_num      -0.9124     0.4935  -1.849  0.06635 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 100.7 on 160 degrees of freedom
Multiple R-squared:  0.3929,    Adjusted R-squared:  0.3815 
F-statistic: 34.52 on 3 and 160 DF,  p-value: < 2.2e-16

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: observed_otus
           Df  Sum Sq Mean Sq F value    Pr(>F)    
Urban       1  534368  534368  28.009 5.994e-07 ***
Gender      1   69292   69292   3.632   0.05922 .  
Age_num     1   51988   51988   2.725   0.10157    
Residuals 113 2155844   19078                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Urban + Gender + Age_num"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-282.28  -99.17  -23.80  102.44  335.26 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 281.1263    25.2675  11.126  < 2e-16 ***
UrbanMedium 135.0262    25.7063   5.253 7.13e-07 ***
GenderM      50.5734    26.3768   1.917   0.0577 .  
Age_num      -1.3304     0.8059  -1.651   0.1016    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 138.1 on 113 degrees of freedom
Multiple R-squared:  0.2332,    Adjusted R-squared:  0.2128 
F-statistic: 11.46 on 3 and 113 DF,  p-value: 1.293e-06

Next we perform the broken stick regression based on the selected final model

stick_models <- list()
for (YY in Yrs) {
    # YY = '2015'
    print(paste0("-------Following is for ", YY, "--------------"))
    Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
        Age_num >= 1)
    fit.seg <- segmented(lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
    summary(fit.seg) %>% print()
    plot(fit.seg$fitted.values, fit.seg$residuals)
    plot(fit.seg$model$observed_otus, fit.seg$fitted.values)
    qqnorm(fit.seg$residuals)
    if (!"segmented" %in% class(fit.seg)) {
        print(paste(YY, BB, AA, "same as lm model", sep = "-"))
    } else {
        anova(fit.seg, lm_models[[YY]]) %>% print()
        davies.test(lm_models[[YY]], seg.Z = ~Age_num, k = 10) %>% print()
        confint.segmented(fit.seg) %>% print()
        stick_models[[YY]] <- fit.seg
        slope(fit.seg) %>% print()
    }
    
    
}
[1] "-------Following is for 2015--------------"

Call:
lm(formula = paste0(AA, " ~ Urban + Gender + Age_num"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-253.91  -66.40   -0.44   55.81  354.88 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 286.3689    16.1831  17.696  < 2e-16 ***
UrbanMedium 149.1367    16.0434   9.296  < 2e-16 ***
GenderM      44.5666    15.8833   2.806  0.00564 ** 
Age_num      -0.9124     0.4935  -1.849  0.06635 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 100.7 on 160 degrees of freedom
Multiple R-squared:  0.3929,    Adjusted R-squared:  0.3815 
F-statistic: 34.52 on 3 and 160 DF,  p-value: < 2.2e-16

[1] "2015-Right_Hand-observed_otus-same as lm model"
[1] "-------Following is for 2016--------------"

Call:
lm(formula = paste0(AA, " ~ Urban + Gender + Age_num"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-282.28  -99.17  -23.80  102.44  335.26 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 281.1263    25.2675  11.126  < 2e-16 ***
UrbanMedium 135.0262    25.7063   5.253 7.13e-07 ***
GenderM      50.5734    26.3768   1.917   0.0577 .  
Age_num      -1.3304     0.8059  -1.651   0.1016    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 138.1 on 113 degrees of freedom
Multiple R-squared:  0.2332,    Adjusted R-squared:  0.2128 
F-statistic: 11.46 on 3 and 113 DF,  p-value: 1.293e-06

[1] "2016-Right_Hand-observed_otus-same as lm model"

data for plots

dat15 <- Work %>% filter(Body_Site == BB, Year == "2015", SampleGroup == "Villagers", 
    Age_num >= 1)
dat16 <- Work %>% filter(Body_Site == BB, Year == "2016", SampleGroup == "Villagers", 
    Age_num >= 1)

dat15 <- cbind(dat15, predict(lm_models[["2015"]], newdata = dat15, interval = "confidence"))
dat16 <- cbind(dat16, predict(lm_models[["2016"]], newdata = dat16, interval = "confidence"))

dat_p[["2015"]][[AA]] = dat15
dat_p[["2016"]][[AA]] = dat16
shannon
BB = "Right_Hand"
AA = "shannon"

Work_dat <- Work %>% filter(Body_Site == BB, SampleGroup == "Villagers", Age_num >= 
    1)
p <- ggplot(Work_dat, aes(x = Age_num, y = get(AA), color = Urban, shape = Gender)) + 
    geom_point(size = 1.5) + scale_shape_manual(values = c(16, 21)) + scale_color_manual(breaks = names(Urban_color), 
    values = Urban_color) + facet_grid(. ~ Year) + theme_jw() + theme(aspect.ratio = 0.8, 
    panel.background = element_rect(color = "black")) + labs(x = "Age (Years)", y = AA, 
    title = BB)
print(p)

# Year 2015
YY = "2015"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
    Age_num >= 1)

# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender + 
    Urban:Gender, data = Work_dat)

# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)

# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num + 
    Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)

anova(regforward.out)
# Year 2016
YY = "2016"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
    Age_num >= 1)

# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender + 
    Urban:Gender, data = Work_dat)

# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)

# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num + 
    Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)

anova(regforward.out)
# Final linear model
lm_models <- list()
for (YY in Yrs) {
    Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
        Age_num >= 1)
    fit.final <- lm(paste0(AA, " ~ Urban + Gender + Age_num"), data = Work_dat)
    print(paste0("--------Below is ANOVA table for year ", YY))
    anova(fit.final) %>% print()
    
    print(paste0("--------Below is coefficient for year ", YY))
    summary(fit.final) %>% print()
    lm_models[[YY]] <- fit.final
}
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: shannon
           Df Sum Sq Mean Sq F value    Pr(>F)    
Urban       1 19.569 19.5685 40.4076 2.058e-09 ***
Gender      1  2.413  2.4132  4.9832   0.02698 *  
Age_num     1  2.366  2.3663  4.8863   0.02849 *  
Residuals 160 77.484  0.4843                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Urban + Gender + Age_num"), data = Work_dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.20037 -0.41410  0.09709  0.48491  1.29455 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  6.287567   0.111834  56.222  < 2e-16 ***
UrbanMedium  0.665392   0.110869   6.002 1.26e-08 ***
GenderM      0.251744   0.109762   2.294   0.0231 *  
Age_num     -0.007539   0.003411  -2.210   0.0285 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6959 on 160 degrees of freedom
Multiple R-squared:  0.2391,    Adjusted R-squared:  0.2248 
F-statistic: 16.76 on 3 and 160 DF,  p-value: 1.621e-09

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: shannon
           Df  Sum Sq Mean Sq F value    Pr(>F)    
Urban       1  13.575 13.5753 11.8432 0.0008121 ***
Gender      1   2.522  2.5217  2.1999 0.1408027    
Age_num     1   7.293  7.2927  6.3622 0.0130500 *  
Residuals 113 129.527  1.1463                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Urban + Gender + Age_num"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.2594 -0.3309  0.1945  0.5960  1.6278 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  6.205022   0.195854  31.682  < 2e-16 ***
UrbanMedium  0.666162   0.199256   3.343  0.00112 ** 
GenderM      0.306865   0.204453   1.501  0.13617    
Age_num     -0.015757   0.006247  -2.522  0.01305 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.071 on 113 degrees of freedom
Multiple R-squared:  0.153, Adjusted R-squared:  0.1305 
F-statistic: 6.802 on 3 and 113 DF,  p-value: 0.000295

Next we perform the broken stick regression based on the selected final model

stick_models <- list()
for (YY in Yrs) {
    # YY = '2015'
    print(paste0("-------Following is for ", YY, "--------------"))
    Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers", 
        Age_num >= 1)
    fit.seg <- segmented(lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
    summary(fit.seg) %>% print()
    plot(fit.seg$fitted.values, fit.seg$residuals)
    plot(fit.seg$model$shannon, fit.seg$fitted.values)
    qqnorm(fit.seg$residuals)
    if (!"segmented" %in% class(fit.seg)) {
        print(paste(YY, BB, AA, "same as lm model", sep = "-"))
    } else {
        anova(fit.seg, lm_models[[YY]]) %>% print()
        davies.test(lm_models[[YY]], seg.Z = ~Age_num, k = 10) %>% print()
        confint.segmented(fit.seg) %>% print()
        stick_models[[YY]] <- fit.seg
        slope(fit.seg) %>% print()
    }
    
    
}
[1] "-------Following is for 2015--------------"

    ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)

Estimated Break-Point(s):
               Est. St.Err
psi1.Age_num 4.314   1.83

Meaningful coefficients of the linear terms:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   5.7584     0.3410  16.889  < 2e-16 ***
UrbanMedium   0.6897     0.1109   6.217 4.32e-09 ***
GenderM       0.2639     0.1092   2.416   0.0168 *  
Age_num       0.1338     0.1195   1.120   0.2643    
U1.Age_num   -0.1445     0.1195  -1.210       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6912 on 158 degrees of freedom
Multiple R-Squared: 0.2587,  Adjusted R-squared: 0.2352 

Convergence attained in 2 iter. (rel. change 1.8799e-16)

Analysis of Variance Table

Model 1: shannon ~ Urban + Gender + Age_num + U1.Age_num + psi1.Age_num
Model 2: shannon ~ Urban + Gender + Age_num
  Res.Df    RSS Df Sum of Sq     F Pr(>F)
1    158 75.492                          
2    160 77.484 -2   -1.9924 2.085 0.1277

    Davies' test for a change in the slope

data:  formula = shannon ~ Urban + Gender + Age_num ,   method = lm 
model = gaussian , link = identity  
segmented variable = Age_num
'best' at = 7.5556, n.points = 9, p-value = 0.4543
alternative hypothesis: two.sided

                Est. CI(95%).low CI(95%).up
psi1.Age_num 4.31412    0.699437     7.9288
$Age_num
            Est.   St.Err. t value CI(95%).l  CI(95%).u
slope1  0.133820 0.1194600  1.1202 -0.102120  0.3697600
slope2 -0.010718 0.0038693 -2.7701 -0.018361 -0.0030762

[1] "-------Following is for 2016--------------"

    ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)

Estimated Break-Point(s):
              Est. St.Err
psi1.Age_num   45  8.018

Meaningful coefficients of the linear terms:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  6.34024    0.21078  30.080  < 2e-16 ***
UrbanMedium  0.69974    0.19811   3.532 0.000602 ***
GenderM      0.27327    0.20319   1.345 0.181383    
Age_num     -0.02709    0.00933  -2.904 0.004452 ** 
U1.Age_num   0.07724    0.04678   1.651       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.06 on 111 degrees of freedom
Multiple R-Squared: 0.184,  Adjusted R-squared: 0.1473 

Convergence attained in 1 iter. (rel. change 2.109e-07)

Analysis of Variance Table

Model 1: shannon ~ Urban + Gender + Age_num + U1.Age_num + psi1.Age_num
Model 2: shannon ~ Urban + Gender + Age_num
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    111 124.78                           
2    113 129.53 -2   -4.7481 2.1119 0.1258

    Davies' test for a change in the slope

data:  formula = shannon ~ Urban + Gender + Age_num ,   method = lm 
model = gaussian , link = identity  
segmented variable = Age_num
'best' at = 46.889, n.points = 9, p-value = 0.2338
alternative hypothesis: two.sided

                Est. CI(95%).low CI(95%).up
psi1.Age_num 44.9996     29.1121    60.8871
$Age_num
            Est.   St.Err. t value CI(95%).l  CI(95%).u
slope1 -0.027091 0.0093302 -2.9036 -0.045580 -0.0086029
slope2  0.050152 0.0458170  1.0946 -0.040638  0.1409400

data for plots

dat15 <- Work %>% filter(Body_Site == BB, Year == "2015", SampleGroup == "Villagers", 
    Age_num >= 1)
dat16 <- Work %>% filter(Body_Site == BB, Year == "2016", SampleGroup == "Villagers", 
    Age_num >= 1)

dat15 <- cbind(dat15, predict(lm_models[["2015"]], newdata = dat15, interval = "confidence"))
dat16 <- cbind(dat16, predict(lm_models[["2016"]], newdata = dat16, interval = "confidence"))

dat_p[["2015"]][[AA]] = dat15
dat_p[["2016"]][[AA]] = dat16
Data for right hand
dat_all[[BB]] = dat_p

Stats check

Stats on the urban, age, and gender

Indices <- colnames(alpha[c(2, 4:5)])
Yrs <- c("2015", "2016")
BodySites <- c("Feces", "Mouth", "Nose", "Right_Arm", "Right_Hand")

lm_results <- list()
duncan_results <- list()

for (AA in Indices) {
    for (BB in BodySites) {
        for (YY in Yrs) {
            Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == 
                "Villagers", Age_num >= 1)
            fit.final <- lm(paste0(AA, " ~ Urban + Gender + Age_num"), data = Work_dat)
            print(paste0("--------Below is ANOVA table for year ", YY))
            anova(fit.final) %>% print()
            
            print(paste0("--------Below is coefficient for year ", YY))
            summary(fit.final) %>% print()
            lm_results[[paste(AA, BB, YY, sep = "-")]] <- summary(fit.final)$coefficients[-1, 
                ] %>% as.data.frame() %>% rownames_to_column(var = "Factor") %>% 
                mutate(Year = YY, Body_Sites = BB, Index = AA)
            fit.final.village <- lm(paste0(AA, " ~ Village + Gender + Age_num"), 
                data = Work_dat)
            duncan_results[[paste(AA, BB, YY, sep = "-")]] <- duncan.test(fit.final.village, 
                trt = "Village", console = T)$groups %>% rownames_to_column(var = "Village") %>% 
                mutate(Year = YY, Body_Sites = BB, Index = AA)
            colnames(duncan_results[[paste(AA, BB, YY, sep = "-")]])[2] <- "Value"
        }
    }
}

lm_results_final <- do.call("rbind", lm_results)
lm_results_final$fdr <- p.adjust(lm_results_final$`Pr(>|t|)`, method = "fdr")
duncan_results_final <- do.call("rbind", duncan_results)

write.table(lm_results_final, file = "../data/alpha_lm_stats.txt", sep = "\t", row.names = F, 
    quote = F)

# clip <- pipe('pbcopy', 'w') write.table(lm_results_final, file = clip,
# row.names = F, quote = F, sep = '\t') close(clip)

within Sanema village effect

Indices <- colnames(alpha[c(2, 4:5)])
Yrs <- c("2015", "2016")
BodySites <- c("Feces", "Mouth", "Nose", "Right_Arm", "Right_Hand")

lm_results <- list()
anova_results <- list()
for (AA in Indices) {
    for (BB in BodySites) {
        for (YY in Yrs) {
            Work_dat <- Work %>% filter(Ethnicity == "SANEMA", Body_Site == BB, Year == 
                YY, SampleGroup == "Villagers", Age_num >= 1)
            fit.final <- lm(paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)
            print(paste0("--------Below is ANOVA table for year ", YY))
            anova(fit.final) %>% print()
            
            print(paste0("--------Below is coefficient for year ", YY))
            summary(fit.final) %>% print()
            lm_results[[paste(AA, BB, YY, sep = "-")]] <- summary(fit.final)$coefficients[-1, 
                ] %>% as.data.frame() %>% rownames_to_column(var = "Factor") %>% 
                mutate(Year = YY, Body_Sites = BB, Index = AA)
            anova_results[[paste(AA, BB, YY, sep = "-")]] <- as.data.frame(anova(fit.final))[1:3, 
                4:5] %>% rownames_to_column(var = "Factor") %>% mutate(Year = YY, 
                Body_Sites = BB, Index = AA)
            
        }
    }
}
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: faith_pd
          Df  Sum Sq Mean Sq F value Pr(>F)  
Village    4  267.74  66.936  1.9358 0.1119  
Gender     1  146.23 146.232  4.2292 0.0428 *
Age_num    1   61.98  61.976  1.7924 0.1842  
Residuals 85 2939.03  34.577                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-14.679  -4.596   1.244   4.749   9.781 

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)    
(Intercept)                30.05499    2.08522  14.413   <2e-16 ***
VillageKuyuwininha         -3.45478    2.24784  -1.537   0.1280    
VillageMosenahanha         -0.04427    2.55933  -0.017   0.9862    
VillageShianana-Jiyakwanha -3.08565    2.18967  -1.409   0.1624    
VillageWashudihanha        -4.49443    2.27934  -1.972   0.0519 .  
GenderM                    -2.90774    1.28950  -2.255   0.0267 *  
Age_num                     0.05008    0.03741   1.339   0.1842    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.88 on 85 degrees of freedom
Multiple R-squared:  0.1394,    Adjusted R-squared:  0.07862 
F-statistic: 2.294 on 6 and 85 DF,  p-value: 0.04222

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: faith_pd
          Df  Sum Sq Mean Sq F value Pr(>F)
Village    4  127.26  31.814  1.0326 0.3964
Gender     1   59.86  59.857  1.9429 0.1676
Age_num    1   17.78  17.783  0.5772 0.4499
Residuals 72 2218.22  30.809               
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-12.5271  -3.8100  -0.1492   3.3117  15.9465 

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)    
(Intercept)                25.03124    2.17450  11.511   <2e-16 ***
VillageKuyuwininha         -2.67620    2.22460  -1.203    0.233    
VillageShianana-Jiyakwanha -2.44842    2.32560  -1.053    0.296    
VillageSudukuma            -1.11259    2.42974  -0.458    0.648    
VillageWashudihanha         1.66737    3.03047   0.550    0.584    
GenderM                     1.66382    1.29812   1.282    0.204    
Age_num                     0.02840    0.03738   0.760    0.450    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.551 on 72 degrees of freedom
Multiple R-squared:  0.08456,   Adjusted R-squared:  0.008272 
F-statistic: 1.108 on 6 and 72 DF,  p-value: 0.366

[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: faith_pd
          Df Sum Sq Mean Sq F value    Pr(>F)    
Village    4  32.32   8.080  0.7172    0.5824    
Gender     1   3.37   3.370  0.2991    0.5858    
Age_num    1 222.48 222.483 19.7495 2.589e-05 ***
Residuals 87 980.08  11.265                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.2432 -1.7891 -0.2504  1.5857 10.1863 

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)    
(Intercept)                12.10518    1.54323   7.844 1.03e-11 ***
VillageKuyuwininha         -1.41335    1.53533  -0.921    0.360    
VillageMosenahanha          0.08616    1.69734   0.051    0.960    
VillageShianana-Jiyakwanha -0.40834    1.53221  -0.267    0.790    
VillageWashudihanha         0.32241    1.59235   0.202    0.840    
GenderM                    -0.10671    0.71749  -0.149    0.882    
Age_num                     0.09539    0.02147   4.444 2.59e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.356 on 87 degrees of freedom
Multiple R-squared:  0.2085,    Adjusted R-squared:  0.1539 
F-statistic:  3.82 on 6 and 87 DF,  p-value: 0.001995

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: faith_pd
          Df Sum Sq Mean Sq F value    Pr(>F)    
Village    3   6.46    2.15  0.1136   0.95172    
Gender     1  60.86   60.86  3.2110   0.07919 .  
Age_num    1 507.04  507.04 26.7513 4.091e-06 ***
Residuals 50 947.69   18.95                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)

Residuals:
   Min     1Q Median     3Q    Max 
-9.407 -3.232 -0.059  2.255 11.458 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)         10.12440    1.53300   6.604 2.50e-08 ***
VillageKuyuwininha  -0.51548    1.61437  -0.319    0.751    
VillageSudukuma      1.18120    1.83064   0.645    0.522    
VillageWashudihanha -0.94686    1.97223  -0.480    0.633    
GenderM              0.06533    1.36342   0.048    0.962    
Age_num              0.19242    0.03720   5.172 4.09e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.354 on 50 degrees of freedom
Multiple R-squared:  0.3774,    Adjusted R-squared:  0.3151 
F-statistic: 6.061 on 5 and 50 DF,  p-value: 0.0001848

[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: faith_pd
          Df Sum Sq Mean Sq F value   Pr(>F)   
Village    4 1439.7  359.92  3.7874 0.007236 **
Gender     1   20.3   20.33  0.2139 0.644995   
Age_num    1  597.2  597.17  6.2839 0.014264 * 
Residuals 78 7412.4   95.03                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-16.839  -6.341  -1.912   4.169  25.992 

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)  
(Intercept)                 7.33590    4.63453   1.583   0.1175  
VillageKuyuwininha         11.12761    4.78519   2.325   0.0226 *
VillageMosenahanha          1.96005    5.29088   0.370   0.7120  
VillageShianana-Jiyakwanha  6.26461    4.82051   1.300   0.1976  
VillageWashudihanha        10.13580    4.94157   2.051   0.0436 *
GenderM                    -2.13749    2.19252  -0.975   0.3326  
Age_num                     0.16052    0.06403   2.507   0.0143 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.748 on 78 degrees of freedom
Multiple R-squared:  0.2172,    Adjusted R-squared:  0.157 
F-statistic: 3.608 on 6 and 78 DF,  p-value: 0.003288

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: faith_pd
          Df Sum Sq Mean Sq F value  Pr(>F)  
Village    3 1000.7  333.57  3.3757 0.02678 *
Gender     1   14.0   13.95  0.1412 0.70895  
Age_num    1  490.8  490.77  4.9665 0.03112 *
Residuals 43 4249.0   98.81                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-17.833  -6.404  -1.736   4.659  23.402 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)          15.9428     3.8103   4.184 0.000139 ***
VillageKuyuwininha    3.0606     4.0422   0.757 0.453084    
VillageSudukuma      -4.9909     4.4679  -1.117 0.270172    
VillageWashudihanha   4.6992     5.7437   0.818 0.417782    
GenderM              -3.6879     3.3868  -1.089 0.282258    
Age_num               0.2237     0.1004   2.229 0.031122 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.941 on 43 degrees of freedom
Multiple R-squared:  0.2616,    Adjusted R-squared:  0.1758 
F-statistic: 3.047 on 5 and 43 DF,  p-value: 0.01928

[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: faith_pd
          Df Sum Sq Mean Sq F value    Pr(>F)    
Village    4 2227.4  556.86  5.6880 0.0004135 ***
Gender     1  388.1  388.08  3.9640 0.0496575 *  
Age_num    1   40.9   40.94  0.4182 0.5195785    
Residuals 86 8419.4   97.90                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-19.9594  -5.8154  -0.0413   4.7356  30.0663 

Coefficients:
                            Estimate Std. Error t value Pr(>|t|)    
(Intercept)                 54.70341    4.14874  13.186  < 2e-16 ***
VillageKuyuwininha         -16.36338    4.24321  -3.856 0.000222 ***
VillageMosenahanha         -17.91609    4.64464  -3.857 0.000221 ***
VillageShianana-Jiyakwanha -15.03259    4.22682  -3.556 0.000614 ***
VillageWashudihanha        -10.90276    4.47693  -2.435 0.016944 *  
GenderM                      3.98777    2.11464   1.886 0.062699 .  
Age_num                      0.03924    0.06069   0.647 0.519579    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.894 on 86 degrees of freedom
Multiple R-squared:  0.2398,    Adjusted R-squared:  0.1868 
F-statistic: 4.522 on 6 and 86 DF,  p-value: 0.0004973

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: faith_pd
          Df Sum Sq Mean Sq F value Pr(>F)  
Village    3 2057.4  685.81  3.8010 0.0220 *
Gender     1   22.3   22.28  0.1235 0.7281  
Age_num    1  104.1  104.10  0.5770 0.4543  
Residuals 26 4691.1  180.43                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-26.215  -8.914   3.249   8.816  31.260 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)          52.7982     6.2340   8.469 5.97e-09 ***
VillageKuyuwininha  -27.0667     9.3237  -2.903  0.00744 ** 
VillageSudukuma      -5.3180     6.1586  -0.864  0.39575    
VillageWashudihanha -17.0004     7.5285  -2.258  0.03256 *  
GenderM              -1.5252     5.9224  -0.258  0.79880    
Age_num               0.1276     0.1680   0.760  0.45434    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 13.43 on 26 degrees of freedom
Multiple R-squared:  0.3176,    Adjusted R-squared:  0.1864 
F-statistic: 2.421 on 5 and 26 DF,  p-value: 0.06271

[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: faith_pd
          Df Sum Sq Mean Sq F value   Pr(>F)   
Village    4  976.1  244.02  3.8654 0.006205 **
Gender     1  435.4  435.36  6.8964 0.010242 * 
Age_num    1    7.7    7.70  0.1220 0.727727   
Residuals 85 5366.0   63.13                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-14.9620  -4.7727  -0.1503   4.2880  24.7994 

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)    
(Intercept)                44.60976    3.49746  12.755  < 2e-16 ***
VillageKuyuwininha         -9.63597    3.64022  -2.647  0.00967 ** 
VillageMosenahanha         -8.83419    3.91219  -2.258  0.02650 *  
VillageShianana-Jiyakwanha -4.03130    3.61637  -1.115  0.26811    
VillageWashudihanha        -6.79309    3.77159  -1.801  0.07523 .  
GenderM                     4.52324    1.70738   2.649  0.00962 ** 
Age_num                    -0.01731    0.04955  -0.349  0.72773    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.945 on 85 degrees of freedom
Multiple R-squared:  0.2092,    Adjusted R-squared:  0.1533 
F-statistic: 3.747 on 6 and 85 DF,  p-value: 0.002348

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: faith_pd
          Df Sum Sq Mean Sq F value   Pr(>F)   
Village    3 2271.4  757.12  5.7523 0.002796 **
Gender     1  101.5  101.49  0.7711 0.386227   
Age_num    1   10.2   10.19  0.0774 0.782559   
Residuals 33 4343.5  131.62                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-23.990  -5.218  -1.464   7.845  20.508 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)          44.59878    3.99891  11.153 9.84e-13 ***
VillageKuyuwininha  -22.14827    6.14464  -3.604  0.00102 ** 
VillageSudukuma       0.34254    4.40263   0.078  0.93845    
VillageWashudihanha  -5.75690    6.34527  -0.907  0.37084    
GenderM               3.26711    4.06429   0.804  0.42723    
Age_num               0.03458    0.12427   0.278  0.78256    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 11.47 on 33 degrees of freedom
Multiple R-squared:  0.3543,    Adjusted R-squared:  0.2564 
F-statistic: 3.621 on 5 and 33 DF,  p-value: 0.01013

[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: shannon
          Df Sum Sq Mean Sq F value Pr(>F)
Village    4  2.993 0.74831  1.0878 0.3678
Gender     1  1.363 1.36264  1.9809 0.1629
Age_num    1  0.450 0.44994  0.6541 0.4209
Residuals 85 58.471 0.68789               
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.0824 -0.6423  0.1359  0.6592  1.3825 

Coefficients:
                            Estimate Std. Error t value Pr(>|t|)    
(Intercept)                 5.870701   0.294116  19.961   <2e-16 ***
VillageKuyuwininha         -0.174672   0.317054  -0.551    0.583    
VillageMosenahanha          0.310353   0.360989   0.860    0.392    
VillageShianana-Jiyakwanha -0.024449   0.308848  -0.079    0.937    
VillageWashudihanha        -0.215965   0.321496  -0.672    0.504    
GenderM                    -0.277349   0.181881  -1.525    0.131    
Age_num                     0.004267   0.005276   0.809    0.421    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.8294 on 85 degrees of freedom
Multiple R-squared:  0.07595,   Adjusted R-squared:  0.01072 
F-statistic: 1.164 on 6 and 85 DF,  p-value: 0.3331

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: shannon
          Df Sum Sq Mean Sq F value Pr(>F)
Village    4  3.922 0.98047  1.9663 0.1088
Gender     1  0.314 0.31397  0.6297 0.4301
Age_num    1  0.000 0.00001  0.0000 0.9968
Residuals 72 35.902 0.49864               
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.84682 -0.42938  0.02862  0.48796  1.43152 

Coefficients:
                             Estimate Std. Error t value Pr(>|t|)    
(Intercept)                 5.460e+00  2.766e-01  19.736   <2e-16 ***
VillageKuyuwininha         -2.237e-01  2.830e-01  -0.791    0.432    
VillageShianana-Jiyakwanha  2.524e-01  2.959e-01   0.853    0.396    
VillageSudukuma            -1.078e-01  3.091e-01  -0.349    0.728    
VillageWashudihanha         3.549e-01  3.855e-01   0.920    0.360    
GenderM                     1.298e-01  1.651e-01   0.786    0.434    
Age_num                     1.933e-05  4.755e-03   0.004    0.997    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7061 on 72 degrees of freedom
Multiple R-squared:  0.1055,    Adjusted R-squared:  0.03099 
F-statistic: 1.416 on 6 and 72 DF,  p-value: 0.2204

[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: shannon
          Df Sum Sq Mean Sq F value Pr(>F)
Village    4  2.976 0.74390  1.1975 0.3177
Gender     1  0.512 0.51212  0.8244 0.3664
Age_num    1  1.091 1.09103  1.7564 0.1885
Residuals 87 54.044 0.62119               
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.35162 -0.36635  0.02832  0.49628  1.73734 

Coefficients:
                            Estimate Std. Error t value Pr(>|t|)    
(Intercept)                 4.293606   0.362386  11.848   <2e-16 ***
VillageKuyuwininha         -0.048434   0.360533  -0.134    0.893    
VillageMosenahanha          0.477654   0.398575   1.198    0.234    
VillageShianana-Jiyakwanha -0.031306   0.359799  -0.087    0.931    
VillageWashudihanha         0.276176   0.373921   0.739    0.462    
GenderM                    -0.185748   0.168483  -1.102    0.273    
Age_num                     0.006680   0.005041   1.325    0.189    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7882 on 87 degrees of freedom
Multiple R-squared:  0.07811,   Adjusted R-squared:  0.01453 
F-statistic: 1.228 on 6 and 87 DF,  p-value: 0.2996

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: shannon
          Df Sum Sq Mean Sq F value    Pr(>F)    
Village    3  2.998  0.9994  1.1125 0.3528618    
Gender     1  1.019  1.0191  1.1344 0.2919507    
Age_num    1 11.602 11.6019 12.9146 0.0007441 ***
Residuals 50 44.917  0.8983                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.97103 -0.65291 -0.08837  0.53085  2.28247 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)          3.387193   0.333748  10.149 9.75e-14 ***
VillageKuyuwininha  -0.271446   0.351462  -0.772 0.443551    
VillageSudukuma      0.675170   0.398546   1.694 0.096471 .  
VillageWashudihanha -0.254288   0.429371  -0.592 0.556362    
GenderM             -0.040747   0.296829  -0.137 0.891366    
Age_num              0.029106   0.008099   3.594 0.000744 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9478 on 50 degrees of freedom
Multiple R-squared:  0.258, Adjusted R-squared:  0.1838 
F-statistic: 3.477 on 5 and 50 DF,  p-value: 0.008957

[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: shannon
          Df Sum Sq Mean Sq F value   Pr(>F)   
Village    4 15.040  3.7601  3.9173 0.005968 **
Gender     1  0.890  0.8903  0.9276 0.338471   
Age_num    1  4.885  4.8855  5.0898 0.026868 * 
Residuals 78 74.869  0.9599                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.6370 -0.5574  0.0151  0.5295  2.3070 

Coefficients:
                            Estimate Std. Error t value Pr(>|t|)    
(Intercept)                 2.556815   0.465776   5.489 4.85e-07 ***
VillageKuyuwininha          0.866755   0.480917   1.802   0.0754 .  
VillageMosenahanha         -0.213381   0.531740  -0.401   0.6893    
VillageShianana-Jiyakwanha  0.318818   0.484468   0.658   0.5124    
VillageWashudihanha         0.709996   0.496634   1.430   0.1568    
GenderM                    -0.311185   0.220351  -1.412   0.1619    
Age_num                     0.014519   0.006436   2.256   0.0269 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9797 on 78 degrees of freedom
Multiple R-squared:  0.2175,    Adjusted R-squared:  0.1574 
F-statistic: 3.614 on 6 and 78 DF,  p-value: 0.003246

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: shannon
          Df Sum Sq Mean Sq F value  Pr(>F)  
Village    3  8.257 2.75239  2.5850 0.06547 .
Gender     1  0.587 0.58698  0.5513 0.46183  
Age_num    1  2.009 2.00935  1.8872 0.17664  
Residuals 43 45.784 1.06474                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.60663 -0.55941  0.05706  0.45811  2.48434 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)          3.21546    0.39552   8.130 3.11e-10 ***
VillageKuyuwininha   0.25954    0.41959   0.619    0.539    
VillageSudukuma     -0.40686    0.46378  -0.877    0.385    
VillageWashudihanha  0.81756    0.59622   1.371    0.177    
GenderM             -0.40555    0.35156  -1.154    0.255    
Age_num              0.01431    0.01042   1.374    0.177    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.032 on 43 degrees of freedom
Multiple R-squared:  0.1916,    Adjusted R-squared:  0.09763 
F-statistic: 2.039 on 5 and 43 DF,  p-value: 0.09217

[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: shannon
          Df Sum Sq Mean Sq F value   Pr(>F)   
Village    4  9.095 2.27378  5.0370 0.001075 **
Gender     1  2.925 2.92506  6.4797 0.012694 * 
Age_num    1  0.015 0.01469  0.0325 0.857273   
Residuals 86 38.822 0.45142                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.23705 -0.31392  0.06807  0.40145  1.70577 

Coefficients:
                             Estimate Std. Error t value Pr(>|t|)    
(Intercept)                 6.9816042  0.2817180  24.782  < 2e-16 ***
VillageKuyuwininha         -0.6076504  0.2881329  -2.109  0.03786 *  
VillageMosenahanha         -1.0233288  0.3153916  -3.245  0.00168 ** 
VillageShianana-Jiyakwanha -0.8307218  0.2870199  -2.894  0.00481 ** 
VillageWashudihanha        -0.2943695  0.3040035  -0.968  0.33561    
GenderM                     0.3587055  0.1435935   2.498  0.01439 *  
Age_num                     0.0007434  0.0041210   0.180  0.85727    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6719 on 86 degrees of freedom
Multiple R-squared:  0.2366,    Adjusted R-squared:  0.1834 
F-statistic: 4.443 on 6 and 86 DF,  p-value: 0.0005815

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: shannon
          Df  Sum Sq Mean Sq F value Pr(>F)
Village    3  4.9432 1.64773  1.8099 0.1701
Gender     1  0.8526 0.85265  0.9366 0.3421
Age_num    1  0.0988 0.09875  0.1085 0.7445
Residuals 26 23.6700 0.91038               
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.9357 -0.3463  0.2187  0.5007  1.2221 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)          6.947236   0.442820  15.689 8.97e-15 ***
VillageKuyuwininha  -1.461267   0.662292  -2.206   0.0364 *  
VillageSudukuma     -0.043154   0.437464  -0.099   0.9222    
VillageWashudihanha -0.357715   0.534773  -0.669   0.5094    
GenderM             -0.420837   0.420687  -1.000   0.3264    
Age_num             -0.003931   0.011935  -0.329   0.7445    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9541 on 26 degrees of freedom
Multiple R-squared:  0.1994,    Adjusted R-squared:  0.04542 
F-statistic: 1.295 on 5 and 26 DF,  p-value: 0.2964

[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: shannon
          Df Sum Sq Mean Sq F value  Pr(>F)  
Village    4  3.490 0.87262  2.2542 0.06992 .
Gender     1  0.759 0.75883  1.9603 0.16512  
Age_num    1  0.288 0.28827  0.7447 0.39059  
Residuals 85 32.904 0.38711                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.76584 -0.36760  0.05787  0.40467  1.23096 

Coefficients:
                            Estimate Std. Error t value Pr(>|t|)    
(Intercept)                 6.356886   0.273876  23.211   <2e-16 ***
VillageKuyuwininha          0.033964   0.285056   0.119   0.9054    
VillageMosenahanha         -0.565078   0.306353  -1.845   0.0686 .  
VillageShianana-Jiyakwanha -0.154536   0.283188  -0.546   0.5867    
VillageWashudihanha        -0.063571   0.295342  -0.215   0.8301    
GenderM                     0.200625   0.133700   1.501   0.1372    
Age_num                    -0.003349   0.003880  -0.863   0.3906    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6222 on 85 degrees of freedom
Multiple R-squared:  0.1212,    Adjusted R-squared:  0.05916 
F-statistic: 1.954 on 6 and 85 DF,  p-value: 0.08139

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: shannon
          Df  Sum Sq Mean Sq F value    Pr(>F)    
Village    3 18.1479  6.0493  7.4850 0.0005919 ***
Gender     1  0.0752  0.0752  0.0931 0.7622392    
Age_num    1  1.7271  1.7271  2.1370 0.1532384    
Residuals 33 26.6704  0.8082                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.68378 -0.64642  0.01611  0.76874  1.56309 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)          6.501204   0.313355  20.747  < 2e-16 ***
VillageKuyuwininha  -1.641425   0.481495  -3.409  0.00174 ** 
VillageSudukuma      0.486595   0.344990   1.410  0.16776    
VillageWashudihanha  0.033095   0.497216   0.067  0.94733    
GenderM             -0.001152   0.318479  -0.004  0.99713    
Age_num             -0.014235   0.009738  -1.462  0.15324    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.899 on 33 degrees of freedom
Multiple R-squared:  0.4279,    Adjusted R-squared:  0.3412 
F-statistic: 4.937 on 5 and 33 DF,  p-value: 0.001753

[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: observed_otus
          Df Sum Sq Mean Sq F value  Pr(>F)  
Village    4  71933 17983.2  2.9929 0.02310 *
Gender     1  23138 23137.8  3.8508 0.05300 .
Age_num    1  16751 16750.6  2.7878 0.09867 .
Residuals 85 510733  6008.6                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-176.646  -55.980    6.271   53.262  178.467 

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)    
(Intercept)                281.7795    27.4882  10.251   <2e-16 ***
VillageKuyuwininha         -51.6439    29.6320  -1.743   0.0850 .  
VillageMosenahanha           7.4068    33.7382   0.220   0.8268    
VillageShianana-Jiyakwanha -27.1231    28.8651  -0.940   0.3501    
VillageWashudihanha        -69.4344    30.0472  -2.311   0.0233 *  
GenderM                    -37.7143    16.9987  -2.219   0.0292 *  
Age_num                      0.8233     0.4931   1.670   0.0987 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 77.52 on 85 degrees of freedom
Multiple R-squared:  0.1796,    Adjusted R-squared:  0.1217 
F-statistic: 3.102 on 6 and 85 DF,  p-value: 0.008525

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: observed_otus
          Df Sum Sq Mean Sq F value Pr(>F)
Village    4  19671  4917.7  1.4301 0.2328
Gender     1   9212  9212.3  2.6790 0.1060
Age_num    1     81    80.8  0.0235 0.8786
Residuals 72 247584  3438.7               
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-121.01  -39.47   -2.29   38.95  196.28 

Coefficients:
                            Estimate Std. Error t value Pr(>|t|)    
(Intercept)                213.77144   22.97306   9.305 5.66e-14 ***
VillageKuyuwininha         -35.84516   23.50232  -1.525    0.132    
VillageShianana-Jiyakwanha -14.89448   24.56937  -0.606    0.546    
VillageSudukuma            -17.21544   25.66953  -0.671    0.505    
VillageWashudihanha         15.91064   32.01615   0.497    0.621    
GenderM                     21.97504   13.71427   1.602    0.113    
Age_num                      0.06054    0.39486   0.153    0.879    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 58.64 on 72 degrees of freedom
Multiple R-squared:  0.1047,    Adjusted R-squared:  0.03013 
F-statistic: 1.404 on 6 and 72 DF,  p-value: 0.225

[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: observed_otus
          Df Sum Sq Mean Sq F value  Pr(>F)  
Village    4   5151  1287.7  1.0305 0.39621  
Gender     1    135   134.6  0.1077 0.74358  
Age_num    1   8307  8306.6  6.6474 0.01161 *
Residuals 87 108714  1249.6                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)

Residuals:
   Min     1Q Median     3Q    Max 
-57.02 -25.67  -2.26  20.37 137.47 

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)    
(Intercept)                 90.6802    16.2534   5.579 2.69e-07 ***
VillageKuyuwininha         -16.0785    16.1702  -0.994   0.3228    
VillageMosenahanha           7.9291    17.8765   0.444   0.6585    
VillageShianana-Jiyakwanha -10.3576    16.1373  -0.642   0.5227    
VillageWashudihanha          3.5212    16.7707   0.210   0.8342    
GenderM                     -5.4706     7.5566  -0.724   0.4710    
Age_num                      0.5829     0.2261   2.578   0.0116 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 35.35 on 87 degrees of freedom
Multiple R-squared:  0.1111,    Adjusted R-squared:  0.04983 
F-statistic: 1.813 on 6 and 87 DF,  p-value: 0.1058

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: observed_otus
          Df Sum Sq Mean Sq F value    Pr(>F)    
Village    3    905   301.5  0.1785 0.9104782    
Gender     1   1151  1151.3  0.6816 0.4129697    
Age_num    1  27029 27028.8 16.0015 0.0002091 ***
Residuals 50  84457  1689.1                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-81.170 -31.132  -1.022  17.458 100.691 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)          68.4666    14.4720   4.731 1.87e-05 ***
VillageKuyuwininha   -3.5870    15.2401  -0.235 0.814885    
VillageSudukuma      19.0999    17.2818   1.105 0.274360    
VillageWashudihanha  -7.2194    18.6184  -0.388 0.699843    
GenderM              -6.3581    12.8711  -0.494 0.623483    
Age_num               1.4049     0.3512   4.000 0.000209 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 41.1 on 50 degrees of freedom
Multiple R-squared:  0.2562,    Adjusted R-squared:  0.1818 
F-statistic: 3.444 on 5 and 50 DF,  p-value: 0.009446

[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: observed_otus
          Df Sum Sq Mean Sq F value   Pr(>F)   
Village    4  75831 18957.7  3.8095 0.007003 **
Gender     1   2331  2330.6  0.4683 0.495789   
Age_num    1  18054 18054.5  3.6280 0.060502 . 
Residuals 78 388162  4976.4                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-111.15  -46.94  -12.22   21.41  216.82 

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)  
(Intercept)                 31.5110    33.5376   0.940   0.3503  
VillageKuyuwininha          81.4254    34.6279   2.351   0.0212 *
VillageMosenahanha          13.8680    38.2873   0.362   0.7182  
VillageShianana-Jiyakwanha  39.0188    34.8835   1.119   0.2668  
VillageWashudihanha         72.1992    35.7595   2.019   0.0469 *
GenderM                    -16.9188    15.8661  -1.066   0.2896  
Age_num                      0.8826     0.4634   1.905   0.0605 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 70.54 on 78 degrees of freedom
Multiple R-squared:  0.1986,    Adjusted R-squared:  0.137 
F-statistic: 3.222 on 6 and 78 DF,  p-value: 0.007017

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: observed_otus
          Df Sum Sq Mean Sq F value  Pr(>F)  
Village    3  49014 16337.9  3.6440 0.01988 *
Gender     1   5726  5725.9  1.2771 0.26471  
Age_num    1   9324  9323.8  2.0796 0.15653  
Residuals 43 192791  4483.5                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-113.62  -38.19  -14.26   30.16  213.71 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)   
(Intercept)          88.7345    25.6661   3.457  0.00124 **
VillageKuyuwininha   27.0840    27.2280   0.995  0.32544   
VillageSudukuma     -34.4691    30.0953  -1.145  0.25841   
VillageWashudihanha  54.5309    38.6894   1.409  0.16589   
GenderM             -35.1792    22.8132  -1.542  0.13039   
Age_num               0.9749     0.6760   1.442  0.15653   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 66.96 on 43 degrees of freedom
Multiple R-squared:  0.2494,    Adjusted R-squared:  0.1621 
F-statistic: 2.858 on 5 and 43 DF,  p-value: 0.02582

[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: observed_otus
          Df Sum Sq Mean Sq F value   Pr(>F)   
Village    4 138346   34587  4.6867 0.001806 **
Gender     1  33281   33281  4.5098 0.036573 * 
Age_num    1   1294    1294  0.1753 0.676453   
Residuals 86 634653    7380                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-193.647  -46.368   -6.585   40.051  299.562 

Coefficients:
                            Estimate Std. Error t value Pr(>|t|)    
(Intercept)                 386.9368    36.0200  10.742  < 2e-16 ***
VillageKuyuwininha          -99.5728    36.8402  -2.703 0.008283 ** 
VillageMosenahanha         -145.8888    40.3254  -3.618 0.000501 ***
VillageShianana-Jiyakwanha  -86.3498    36.6979  -2.353 0.020905 *  
VillageWashudihanha         -63.6089    38.8694  -1.636 0.105394    
GenderM                      37.5990    18.3596   2.048 0.043618 *  
Age_num                       0.2206     0.5269   0.419 0.676453    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 85.91 on 86 degrees of freedom
Multiple R-squared:  0.2141,    Adjusted R-squared:  0.1593 
F-statistic: 3.905 on 6 and 86 DF,  p-value: 0.001695

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: observed_otus
          Df Sum Sq Mean Sq F value Pr(>F)  
Village    3 203674   67891  4.3376 0.0132 *
Gender     1   4310    4310  0.2754 0.6042  
Age_num    1  12620   12620  0.8063 0.3775  
Residuals 26 406946   15652                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-222.88  -91.70   27.60   57.89  351.28 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)          420.128     58.063   7.236  1.1e-07 ***
VillageKuyuwininha  -282.647     86.840  -3.255  0.00314 ** 
VillageSudukuma      -43.697     57.360  -0.762  0.45304    
VillageWashudihanha -154.029     70.120  -2.197  0.03716 *  
GenderM              -22.778     55.160  -0.413  0.68303    
Age_num                1.405      1.565   0.898  0.37746    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 125.1 on 26 degrees of freedom
Multiple R-squared:  0.3515,    Adjusted R-squared:  0.2268 
F-statistic: 2.819 on 5 and 26 DF,  p-value: 0.03655

[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table

Response: observed_otus
          Df Sum Sq Mean Sq F value  Pr(>F)  
Village    4  51640 12910.1  2.9234 0.02565 *
Gender     1  15757 15757.1  3.5680 0.06231 .
Age_num    1    120   120.0  0.0272 0.86944  
Residuals 85 375377  4416.2                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"

Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-127.519  -49.928   -0.568   39.532  214.100 

Coefficients:
                            Estimate Std. Error t value Pr(>|t|)    
(Intercept)                305.41513   29.25251  10.441   <2e-16 ***
VillageKuyuwininha         -38.85588   30.44654  -1.276   0.2054    
VillageMosenahanha         -66.76567   32.72127  -2.040   0.0444 *  
VillageShianana-Jiyakwanha  -4.86418   30.24703  -0.161   0.8726    
VillageWashudihanha        -29.78978   31.54525  -0.944   0.3477    
GenderM                     27.05143   14.28041   1.894   0.0616 .  
Age_num                     -0.06833    0.41447  -0.165   0.8694    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 66.45 on 85 degrees of freedom
Multiple R-squared:  0.1524,    Adjusted R-squared:  0.09262 
F-statistic: 2.548 on 6 and 85 DF,  p-value: 0.02564

[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table

Response: observed_otus
          Df Sum Sq Mean Sq F value   Pr(>F)   
Village    3 233404   77801  6.4756 0.001438 **
Gender     1   3103    3103  0.2582 0.614710   
Age_num    1    365     365  0.0304 0.862612   
Residuals 33 396478   12014                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"

Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-219.109  -56.237   -0.805   67.591  248.268 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)          348.3303    38.2059   9.117 1.55e-10 ***
VillageKuyuwininha  -207.8651    58.7065  -3.541  0.00121 ** 
VillageSudukuma       32.8143    42.0632   0.780  0.44088    
VillageWashudihanha  -41.4571    60.6233  -0.684  0.49885    
GenderM               20.6938    38.8307   0.533  0.59766    
Age_num               -0.2071     1.1873  -0.174  0.86261    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 109.6 on 33 degrees of freedom
Multiple R-squared:  0.374, Adjusted R-squared:  0.2792 
F-statistic: 3.943 on 5 and 33 DF,  p-value: 0.006508
lm_results_final <- do.call("rbind", lm_results)
anova_results_final <- do.call("rbind", anova_results)

# clip <- pipe('pbcopy', 'w') write.table(anova_results_final, file = clip,
# row.names = F, quote = F, sep = '\t') close(clip)